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A  GENERALIZED  FINITE -DIFFERENCE  FORMULATION  FOR  THE 
U.S.  GEOLOGICAL  SURVEY  MODULAR  THREE-DIMENSIONAL 
FINITE-DIFFERENCE  GROUND-WATER  FLOW  MODEL 


By  Arlen  W.  Harbaugh 
ABSTRACT 

The  U.S.  Geological  Survey’s  Modular  Ground- Water  Flow  Model  assumes  that  model  nodes  are 
in  the  center  of  cells  and  that  transmissivity  is  constant  within  a  cell.  Based  on  these  assumptions, 
the  model  calculates  coefficients,  called  conductance,  that  are  multiplied  by  head  difference  to 
determine  flow  between  cells.  Although  these  are  common  assumptions  in  finite-difference 
models,  other  assumptions  are  possible.  A  new  option  to  the  model  program  reads  conductance  as 
input  data  rather  than  calculating  it.  This  option  allows  the  user  to  calculate  conductance  outside 
of  the  model.  The  user  thus  has  the  flexibility  to  define  conductance  using  any  desired 
assumptions. 

For  a  water-table  condition,  horizontal  conductance  must  change  as  water  level  varies.  To  handle 
this  situation,  the  new  option  reads  conductance  divided  by  thickness  (CDT)  as  input  data.  The 
model  calculates  saturated  thickness  and  multiplies  it  by  CDT  to  obtain  conductance.  Thus,  the 
user  is  still  free  from  the  assumptions  of  centered  nodes  and  constant  transmissivity  in  cells. 

The  model  option  is  written  in  FORTRAN  77  and  is  fully  compatible  with  the  existing  model.  This 
report  documents  the  new  model  option;  it  includes  a  description  of  the  concepts,  detailed  input 
instructions,  and  a  listing  of  the  code. 


INTRODUCTION 

The  U.S.  Geological  Survey  developed  a  modular  ground-water  flow  simulation  program  that  uses 
the  finite-difference  approximation  (McDonald  and  Harbaugh,  1988).  Although  the  model  can 
simulate  a  variety  of  problems,  it  also  was  designed  so  that  new  capabilities  could  be  easily  added 
if  needed.  This  report  describes  a  new  option  to  the  program  that  removes  some  of  the  assumptions 
that  were  made  for  the  calculation  of  flow  between  model  nodes. 

The  modular  structure  of  the  model  makes  it  easy  to  modify  (McDonald  and  Harbaugh,  1988, 
p.  1-2).  Each  module  is  small  and  compact,  and  the  modules  are  grouped,  mainly  by  hydrologic 
function,  into  packages.  Modules  dealing  with  internal  aquifer  flow,  which  is  flow  between  cells 
and  flow  into  storage,  are  grouped  into  a  package,  herein  called  the  internal  flow  package.  There 
can  be  only  one  internal  flow  package  in  use  at  a  time,  and  in  fact  there  is  only  one  such  package 
in  the  original  model,  the  Block-Centered  Flow  (BCF)  Package. 

To  implement  the  new  option  described  in  this  report,  a  second  internal  flow  package  was 
developed.  The  new  package  is  called  the  General  Finite-Difference  Flow  (GFD)  Package.  The 
structure  of  the  GFD  Package  is  consistent  with  that  of  the  original  model  and  is  written  in  the 
FORTRAN  77  programming  language  (ANSI,  1978)  as  is  the  original  model.  The  BCF  Package 
was  left  in  the  model  so  that  it  may  still  be  used  when  desired.  All  other  model  packages  are 
unaffected  by  the  new  GFD  Package. 
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The  remainder  of  this  report  describes  how  the  GFD  Package  is  implemented  and  how  to  use  it; 
however,  it  is  not  the  purpose  of  this  report  to  describe  when  to  use  the  new  package.  Readers 
should  use  hydrologic  judgement  to  decide  what  assumptions  are  best  for  their  applications.  This 
report  is  intended  to  be  used  as  a  supplement  to  the  documentation  of  the  model  by  McDonald  and 
Harbaugh  (1988). 


The  finite-difference  flow  equation  for  a  model  node  (i,j,k),  where  i,  j,  and  k  are  row,  column,  and 
layer  grid  indexes  respectively,  is  (McDonald  and  Harbaugh,  1988,  p.  2-26): 
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Equation  1  is  written  for  each  time  step  in  a  simulation,  and  the  m  superscript  designates  the  time 
step  number.  Hydraulic  head  is  designated  by  h.  The  HCOF  and  RHS  terms  incorporate  storage 
and  external  stresses.  The  CR,  CC,  and  CV  terms  are  conductances.  Conductance  is  the  product 
of  hydraulic  conductivity  and  cross-sectional  area  of  flow  divided  by  the  length  of  the  flow  path 
and  results  from  application  of  Darcy’s  law  describing  one-dimensional,  steady-state  flow 
occurring  in  a  volume  of  porous  material.  Darcy’s  law  says  that  flow  through  a  volume  of  porous 
material  is  the  product  of  conductance  and  the  head  difference  along  the  flow  path.  The  terms  in 
equation  1  involving  conductance  specify  the  flow  between  the  node  for  which  the  equation  is 
being  written  and  its  six  neighboring  nodes. 


Conductance  of  a  volume  of  porous  material  is  defined  as  (fig.  1) 
C  =  KWb/L 


(2) 


where 


K  is  hydraulic  conductivity  in  the  direction  of  flow, 

W  is  the  width  of  the  volume  in  one  direction  perpendicular  to  flow, 
b  is  the  width  of  the  volume  in  the  other  direction  perpendicular  to  flow,  and 
L  is  the  length  of  the  volume  in  the  direction  of  flow. 

Note  that  the  direction  of  flow  plays  an  important  part  in  the  definition  of  conductance.  The 
conductance  for  a  volume  of  material  may  vary  significantly  for  different  directions  of  flow. 

Each  conductance  in  equation  1  represents  the  properties  of  flow  of  the  material  between  two 
nodes.  That  is,  model  conductance  does  not  apply  to  a  single  cell.  Subscript  indexes  for 
conductance  values  contain  +  or  -1/2  to  indicate  the  area  represented  by  the  conductance.  For 
example,  CC(i-l/2,j,k)  represents  conductance  between  nodes  (i-l,j,k)  and  (i,j,k).  Similarly, 
CR(i,j+l/2,k)  is  the  conductance  between  nodes  (i,j,k)  and  (i,j+l,k).  Figure  2  shows  two  model 
cells,  (i,j,k)  and  (i,j+l,k),  and  the  volume  of  porous  material  represented  by  conductance 
CR(i,j+l/2,k).  Similar  diagrams  could  be  drawn  to  illustrate  model  conductance  between  cell 
(i,j,k)  and  the  other  five  adjacent  cells. 
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hi  hi 


KWb 

C  =  — ; —  where  K  is  hydraulic  conductivity. 


Darcy’s  Law:  Q  =  C  (  h2  -  hi )  where  Q  is  the  flow  through  the  block  and 

h2  and  hi  are  heads  at  the  ends  of  the  block. 

Figure  1.  --  Definition  of  conductance  (C)  for  a  volume  of  porous  material. 


Plan  View 


Explanation 

Node  within  cell. 

Volume  of  material 
represented  by  conductance 
CR(i,j+1/2,k). 


Figure  2.  --  Conductance  between  two  model  nodes. 
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The  BCF  Package  calculates  horizontal  conductance  from  the  hydraulic  conductivity  and  cell 
dimensions  specified  for  individual  cells.  In  order  to  develop  the  equations  for  calculating 
conductances,  it  is  assumed  that  model  nodes  are  located  at  the  center  of  cells,  called  a  block- 
centered  grid,  and  that  transmissivity  is  constant  within  a  cell  (McDonald  and  Harbaugh,  1988, 
p.  5-6).  Although  transmissivity  is  assumed  to  be  constant  within  a  cell,  it  may  vary  from  cell  to 
cell,  which  results  in  discrete  changes  in  transmissivity  at  cell  boundaries. 

The  assumptions  of  a  block-centered  grid  and  constant  transmissivity  within  a  cell  are  not  required 
for  the  development  of  equation  1,  and  the  GFD  Package  removes  these  assumptions.  The  goal  for 
this  new  option  was  to  allow  the  user  full  freedom  to  specify  conductances  as  desired.  For  confined 
cells  this  goal  is  easily  met;  conductance  is  simply  read  as  input  data.  Thus  the  user  now  can 
calculate  conductance  externally  using  any  desired  assumptions.  For  unconfined  cells,  however, 
this  goal  is  only  partly  achievable. 

For  unconfined  cells,  conductance  cannot  be  read  as  input  because  it  must  change  during  the 
simulation  according  to  changes  in  calculated  water  level.  The  GFD  Package  calculates  saturated 
thickness  and  reads  as  input  conductance  divided  by  thickness  between  nodes.  The  program  then 
calculates  conductance  as  the  product  of  saturated  thickness  and  conductance  divided  by  thickness. 
Conductance  divided  by  thickness  represents  material  between  nodes  and  can  be  calculated  in  any 
way  the  user  prefers,  which  is  less  restrictive  than  the  original  model. 

The  GFD  Package  can  be  useful  in  many  situations,  a  few  of  which  are  briefly  mentioned  here. 

1.  A  point-centered  grid  in  which  cell  boundaries  are  equally  spaced  between  nodes  (McDonald 
and  Harbaugh,  1988,  p.  2-5)  is  sometimes  desired  instead  of  a  block-centered  grid. 

2.  The  radial-flow  finite-difference  equation  has  the  same  form  as  equation  1  (Stallman,  1963). 
Thus,  a  radial  flow  model  can  be  constructed  using  the  GFD  Package  if  the  conductances  are 
calculated  using  radial  geometry. 

3.  When  changes  in  transmissivity  between  nodes  are  large,  conductance  values  can  vary 
significantly  depending  on  the  assumptions  regarding  the  distribution  of  transmissivity.  Appel 
(1976)  derives  a  formula  for  calculating  interblock  transmissivity  when  transmissivity  varies 
linearly  between  model  nodes. 

4.  The  geometry  of  an  aquifer  can  sometimes  be  more  accurately  represented  using  the  GFD 
Package.  In  the  BCF  Package,  changing  the  transmissivity  for  one  cell  affects  the  conductances  to 
all  four  neighboring  horizontal  cells.  With  the  GFD  Package,  any  individual  conductance  to  a  cell 
can  be  adjusted  without  impacting  the  other  conductances  to  that  cell.  As  an  example,  this  can  be 
useful  for  simulating  a  flow  barrier;  those  conductances  that  are  affected  by  the  barrier  can  be 
adjusted  without  affecting  other  nearby  conductances. 

5.  Finally,  horizontal  anisotropy  must  be  constant  for  an  entire  layer  in  the  BCF  Package,  but 
the  GFD  Package  can  be  used  to  specify  anisotropy  variations  cell  by  cell.  Note  that  the  principal 
axes  of  hydraulic  conductivity  must  be  aligned  with  the  model  grid,  which  is  an  assumption 
inherent  in  the  partial-differential  flow  equation  on  which  the  model  is  based  (McDonald  and 
Harbaugh,  1988,  p.  2-1). 
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IMPLEMENTATION  OF  THE  GENERAL  FINITE-DIFFERENCE  FLOW  PACKAGE 


The  model  formulates  equation  1  at  each  model  node.  An  internal  flow  package  calculates  the  six 
conductances  (CR,  CC,  and  CV)  used  in  equation  1.  Storage  terms,  which  are  incorporated  in  RHS 
and  HCOF  terms,  are  also  calculated  by  an  internal  flow  package. 

Conductance 

Horizontal  Conductance  for  Confined  Model  Layers 

For  confined  model  cells,  the  horizontal  conductances  (CR  and  CC)  are  read  as  input  data,  and 
therefore  the  user  must  calculate  the  conductance  values  prior  to  making  a  simulation.  These 
conductances  remain  constant  throughout  the  simulation. 

Equation  2  is  used  to  calculate  conductance  for  a  rectangular  grid.  When  calculating  horizontal 
conductance  between  two  nodes  in  a  row  with  a  rectangular  grid,  W  is  the  width  of  the  row,  b  is 
the  thickness  of  the  model  layer  between  the  two  nodes,  L  is  the  horizontal  distance  between  the 
two  nodes,  and  K  is  the  horizontal  hydraulic  conductivity  between  the  nodes.  When  calculating 
horizontal  conductance  between  two  nodes  in  a  column  of  a  rectangular  grid,  W  is  the  width  of  the 
column,  b  is  the  thickness  of  the  model  layer  between  the  two  nodes,  L  is  the  horizontal  distance 
between  the  two  nodes,  and  K  is  the  hydraulic  conductivity  between  the  nodes. 

Horizontal  Conductance  for  Unconfined  and  Convertible  Model  Layers 

In  an  unconfined  aquifer,  horizontal  conductance  changes  as  the  saturated  thickness  of  the  aquifer 
changes.  Thus,  horizontal  conductance  values  must  be  calculated  within  the  model  during  the 
simulation  rather  than  being  specified  as  constants,  as  in  the  confined  situation.  In  the  GFD 
Package,  this  is  accomplished  by  having  the  user  supply  values  of  conductance  divided  by 
thickness  (CDT)  between  nodes.  CDT  is  obtained  by  dividing  equation  2  by  thickness  b: 

CDT  =  KW/L.  (3) 

CDT  values  are  read  as  input  and  represent  the  aquifer  material  between  nodes  just  as  conductance 
values  do.  Given  a  uniform  saturated  thickness  within  a  volume  of  porous  material,  conductance 
is  simply  the  product  of  CDT  and  saturated  thickness.  Saturated  thickness  at  nodes  is  calculated 
within  the  GFD  Package  during  the  simulation  just  as  it  is  done  within  the  BCF  Package 
(McDonald  and  Harbaugh,  1988,  p.  5-9): 


h-BOT 

for 

TOP  >  h  >  BOT, 

(4a) 

TOP-BOT 

for 

h  >  TOP,  or 

(4b) 

0 

for 

h  <  BOT, 

(4c) 

where 

h  is  the  current  head  at  node  (i,j,k), 

BOT  is  the  elevation  of  the  aquifer  bottom  at  node  (i,j,k),  and 
TOP  is  the  elevation  of  the  aquifer  top  at  node  (i,j,k). 
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Saturated  thickness  is  generally  not  uniform  between  nodes,  so  a  formulation  for  calculating 
horizontal  conductance  from  CDT  between  nodes  and  saturated  thickness  at  nodes  must  be 
developed.  The  approach  used  for  the  GFD  Package  is  to  calculate  an  equivalent  saturated 
thickness  that,  when  multiplied  by  CDT  and  head  difference,  results  in  the  correct  flow. 

Appel  (1976)  shows  that  an  equivalent  transmissivity,  T,  for  one- dimensional  flow  in  a  rectangular 
volume  of  porous  material  in  which  transmissivity  varies  linearly  along  the  flow  path  is 

T  =  (T2  -Tl)/ln(T2/Tl)  (5) 


where 


T2  is  transmissivity  at  one  end  of  the  volume  and 
T1  is  transmissivity  at  the  other  end  of  the  volume. 

In  order  to  implement  a  form  of  equation  5  in  the  GFD  Package,  it  is  assumed  that  hydraulic 
conductivity  is  uniform  between  nodes  and  that  the  saturated  thickness  varies  linearly  between 
nodes.  The  product  of  a  constant  hydraulic  conductivity  and  a  linearly  varying  saturated  thickness 
will  produce  a  linearly  varying  transmissivity,  which  is  the  assumption  on  which  equation  5  is 
based.  Clearly,  a  linear  variation  of  saturated  thickness  between  nodes  is  only  an  approximation 
because  of  the  water-table  curvature  and  the  typical  aquifer  bottom  irregularities.  However,  it  is 
generally  possible  to  make  cells  small  enough  (through  specification  of  the  horizontal  grid  spacing) 
so  that  such  an  approximation  is  reasonable. 

If  saturated  thickness  times  a  constant  hydraulic  conductivity  is  substituted  for  transmissivity  in 
equation  (5),  hydraulic  conductivity  cancels,  and  the  result  is 

B  =  (B2-Bl)/ln(B2/Bl)  (6) 


where 


B  is  equivalent  saturated  thickness  for  the  volume, 

B2  is  saturated  thickness  at  one  end  of  the  volume,  and 
B 1  is  saturated  thickness  at  the  other  end  of  the  volume. 

Although  equation  6  provides  a  workable  formulation  of  equivalent  saturated  thickness, 
computational  efficiency  of  the  conductance  calculations  is  also  important  in  the  model  program 
when  there  is  a  large  number  of  unconfined  cells.  Logarithm  calculation  is  computationally 
intensive,  and  an  alternative  is  desirable.  Appel  (1976)  shows  that  when  T1  and  T2  are  nearly 
equal,  the  average  of  T1  and  T2  is  a  good  approximation  of  T  in  equation  5.  It  follows  that 

Bavg  =  (B 1  +  B2)/2,  (7) 

is  a  good  approximation  of  B  in  equation  6  when  B 1  and  B2  are  nearly  equal.  In  particular,  Bavg 
is  within  0.5  percent  of  B  whenever  B2/B 1  is  less  than  1.25  and  greater  than  0.8.  Therefore,  in  the 
GFD  Package,  Bavg  is  used  as  equivalent  saturated  thickness  when  B2/B1  is  less  than  1.25  and 
greater  than  0.8.  Otherwise,  B  as  defined  in  equation  6  is  used. 

The  GFD  Package  calculates  conductance  for  unconfined  layers,  and  layers  that  can  convert 
between  confined  and  unconfined  at  the  beginning  of  each  iteration  using  the  head  from  the 
previous  iteration  to  calculate  saturated  thickness.  If  either  cell  in  this  calculation  has  0  saturated 
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thickness,  then  the  conductance  is  defined  as  0  and  the  cell  becomes  inactive  (is  converted  to  no 
flow).  This  means  that  all  4  horizontal  conductances  to  a  cell  will  be  0  whenever  saturated 
thickness  at  the  cell  is  0.  In  addition,  the  two  vertical  conductances  connecting  to  a  cell  with 
0  saturated  thickness  are  set  to  0.  A  cell  cannot  be  changed  from  inactive  to  active. 

One  further  note  on  the  applicability  of  the  water-table  conductance  calculation  is  that  Appel’s 
(1976)  formulation  for  transmissivity  was  developed  for  a  rectangular  volume  assuming  one¬ 
dimensional  steady-state  flow.  Under  the  assumptions  of  a  water  table  and  a  non-flat  bottom,  the 
volume  being  simulated  by  the  GFD  Package  will  not  be  rectangular  and  flow  will  not  be  one 
dimensional.  Yet,  the  assumptions  are  reasonable  approximations  in  many  situations.  In 
particular,  the  assumptions  are  valid  when  the  aquifer  bottom  variations  and  water-table  curvature 
are  small  compared  to  the  horizontal  grid  spacing.  These  assumptions  are  similar  to,  and  no  more 
restrictive  than,  the  assumptions  made  in  the  BCF  Package  for  calculating  conductance  for 
unconfined  cells;  further,  there  is  no  need  to  make  an  assumption  about  node  location  within  a  cell 
as  in  the  BCF  Package.  Thus  the  goal  of  making  the  GFD  Package  independent  of  node  orientation 
is  maintained. 

The  model  user  must  calculate  CDT  values  outside  of  the  model.  Use  equation  3  to  calculate  CDT 
values  for  a  rectangular  grid.  When  calculating  CDT  between  two  nodes  in  a  row  in  a  rectangular 
grid,  W  is  the  width  of  the  row,  L  is  the  horizontal  distance  between  the  two  nodes,  and  K  is  the 
hydraulic  conductivity  between  the  nodes.  When  calculating  CDT  between  two  nodes  in  a  column 
of  a  rectangular  grid,  W  is  the  width  of  the  column,  L  is  the  horizontal  distance  between  the  two 
nodes,  and  K  is  the  hydraulic  conductivity  between  the  nodes. 

Vertical  Conductance 

Values  of  vertical  conductance  between  cells  (CV)  are  read  as  input  in  all  situations,  whether  the 
connected  nodes  are  specified  as  confined  or  unconfmed.  (Note  that  the  BCF  Package  reads  as 
input  vertical  conductance  divided  by  cell  area  [McDonald  and  Harbaugh,  1988,  p.  5-12]).  In 
unconfined  conditions,  vertical  conductance  does  not  change  with  saturated  thickness;  except  that 
if  a  node  becomes  inactive,  then  both  vertical  conductances  connecting  to  neighboring  nodes  are 
set  to  0.  Use  equation  2  to  calculate  vertical  conductance  values  for  a  rectangular  grid.  In  this 
situation,  W  and  b  are  the  horizontal  dimensions  (input  parameters  DELR  and  DELC)  of  the  two 
cells  containing  the  nodes  between  which  conductance  is  being  calculated,  L  is  the  vertical  distance 
between  the  two  nodes,  and  K  is  the  hydraulic  conductivity  between  the  nodes. 

Vertical  Flow  Calculation  Under  Dewatered  Conditions 


When  a  cell  below  another  cell  becomes  unconfined,  there  is  an  unsaturated  zone  between  them. 
Flow  through  the  unsaturated  zone  is  no  longer  dependent  on  the  head  in  the  lower  cell.  Vertical 
flow  is  a  function  of  the  difference  between  the  head  in  the  upper  cell  and  the  elevation  of  the  top 
of  the  lower  cell.  This  situation  causes  a  change  in  the  basic  finite-difference  equation  1,  which  is 
handled  as  described  by  McDonald  and  Harbaugh  (1988,  p.  5-19).  Compensating  terms  for  the 
erroneous  head  difference  are  added  to  RHS  and  HCOF  coefficients. 

Storage  Terms 

Storage  terms  are  added  to  the  flow  equation  exactly  as  described  by  McDonald  and  Harbaugh 
(1988,  p.  5-24)  for  the  BCF  Package.  A  summary  of  these  terms  follows. 
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For  the  case  of  a  model  cell  that  is  either  strictly  confined  or  strictly  unconfined,  the  expression  for 
the  rate  of  accumulation,  Q,  of  water  in  a  cell  as  a  result  of  storage  is 

Q  =  SCl(h  -  HOLD)/DELT,  (8) 


where 


SCI  is  storage  capacity  (product  of  cell  area  and  specific  yield  or  confined  storage 
coefficient), 

HOLD  is  the  head  at  the  start  of  the  time  step, 
h  is  the  head  at  the  end  of  the  time  step,  and 
DELT  is  the  length  of  the  time  step. 

In  the  formulation  of  equation  (1),  the  term  -SC1/DELT  is  added  to  HCOF  and 
-SC  1  *HOLD/DELT  is  added  to  RHS. 

In  the  case  of  a  cell  that  is  converting  between  confined  and  unconfined  and  vice  versa  during  a 
time  step,  the  expression  for  the  rate  of  accumulation,  Q,  of  water  in  a  cell  as  a  result  of  storage  is 

Q  =  [  SCB(h-TOP)  +  SCA(TOP-HOLD)  ]  /  DELT  (9) 


where 


TOP  is  the  elevation  of  the  top  of  the  aquifer  in  the  cell, 

SCA  is  storage  capacity  in  effect  at  the  start  of  the  time  step,  and 
SCB  is  the  storage  capacity  in  effect  at  the  current  iteration. 

The  term  -SCB/DELT  is  added  to  HCOF  and  term  (-SCB(TOP)+SCA(TOP-HOLD))/DELT  is 
added  to  RHS  in  equation  1. 

Storage  capacity  is  required  for  the  storage  calculations  as  indicated  by  equations  8  and  9.  To  be 
consistent  with  the  goal  of  being  independent  of  the  geometry  of  cells,  storage  capacity  is  directly 
read  as  input  in  the  GFD  Package.  This  is  different  than  in  the  BCF  Package,  where  storage 
coefficient  or  specific  yield  are  read  as  input,  and  the  package  multiples  these  by  cell  area  to  obtain 
storage  capacity. 


APPLICABILITY  AND  LIMITATIONS  OF  OPTIONAL  FORMULATIONS 

The  optional  formulations  for  water-table  conductance,  vertical  leakage  correction,  and  confined/ 
unconfined  storage  conversion  were  developed  using  the  conceptualization  of  a  layered  aquifer 
system  in  which  each  aquifer  is  simulated  by  one  model  layer  and  the  aquifer  layers  are  separated 
by  distinct  confining  units.  Using  these  options  where  these  conditions  are  not  satisfied  may  lead 
to  problems  and  inaccuracies.  For  example,  vertical  conductance  is  left  constant  until  a  cell 
converts  to  no  flow,  and  then  it  is  set  to  zero.  For  this  to  be  reasonable,  there  must  generally  be  a 
confining  layer,  which  is  the  dominant  factor  in  controlling  vertical  flow,  below  the  model  water- 
table  layer.  These  limitations  are  similar  to  the  limitations  for  the  BCF  Package  (McDonald  and 
Harbaugh,  1988,  p.  5-30). 
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DATA  REQUIREMENTS 


The  fundamental  parameters  controlling  cell-to-cell  flow  and  storage  are  read  by  the  General 
Finite-Difference  Package  as  input  data.  These  parameters  may  include  conductance,  conductance 
divided  by  saturated  thickness  (CDT),  storage  capacity,  aquifer  bottom  elevation,  and  aquifer  top 
elevation.  The  specific  parameters  that  are  read  depend  on  which  options  are  invoked. 

The  GFD  Package  is  required  to  define  the  grid  spacing  even  though  the  GFD  Package  does  not 
require  cell  dimensions  in  order  to  calculate  conductance,  storage  capacity,  or  any  other  parameter. 
Although  the  grid  dimensions  are  already  incorporated  into  the  input  parameters  used  by  the  GFD 
Package,  an  internal  flow  package  is  required  to  define  grid  spacing  so  that  it  can  be  used  by  other 
packages.  Other  packages,  for  example  the  recharge  package,  that  use  values  of  grid  spacing  use 
them  only  to  calculate  cell  area.  No  assumptions  about  the  location  of  nodes  within  cells  are 
implied  by  specifying  grid  spacing.  The  grid  spacing  indicates  the  length  of  cells  along  rows 
(DELR)  and  columns  (DELC).  Cell  area  is  calculated  as  the  product  of  DELR  and  DELC  for  that 
cell.  If  non-rectangular  cells  are  being  used,  such  as  would  be  the  case  for  a  radial  flow  model, 
then  values  for  DELR  and  DELC  must  be  specified  so  that  their  product  results  in  the  correct 
cell  area. 

The  user  specifies  a  “layer-type  code”  for  each  layer.  All  cells  in  a  layer  have  the  same  type  code. 
The  layer-type  code  indicates  the  options  in  effect  for  the  cells  in  a  layer,  and  thus  the  data 
requirements  for  a  layer.  Data  for  a  layer  are  input  as  a  series  of  two-dimensional  arrays,  one  array 
for  each  required  parameter.  All  the  required  arrays  are  read  for  a  given  layer  before  reading  any 
data  for  the  next  layer.  That  is,  all  required  arrays  are  read  for  layer  1 ,  then  layer  2  arrays  are  read, 
and  so  on  until  data  has  been  read  for  all  layers. 

The  required  parameters  should  be  specified  for  every  cell  of  a  layer,  including  constant-head  and 
no-flow  cells.  For  no-flow  cells,  the  entered  values  are  never  used  in  calculations;  therefore,  any 
values  may  be  specified.  For  constant-head  cells,  the  storage  capacity  values  are  not  used,  but  the 
other  parameters  are.  Realistic  values  must  be  specified  at  constant-head  cells  for  all  parameters 
except  storage  capacity. 

For  each  node  (i,j,k),  three  conductance  values,  or  two  conductance  divided  by  thickness  (CDT) 
values  and  one  vertical  conductance  value,  are  read  as  input  data.  As  with  all  layer  data, 
conductance  and  CDT  are  read  as  two-dimensional  arrays.  The  names  of  the  conductance  arrays 
correspond  to  the  names  of  variables  in  equation  (1),  which  are  CR,  CC,  and  CV.  The  two  CDT 
arrays  are  named  CDTR  and  CDTC.  Conductances  and  CDT  values  are  indexed  according  to  node 
indexes  even  though  they  represent  material  between  nodes.  The  conductance  between  node  (i,j,k) 
and  node  (i,j+l,k)  is  read  as  CR(i,j,k);  the  conductance  between  node  (i,j,k)  and  node  (i+l,j,k)  is 
read  as  CC(i,j,k).  Likewise,  the  conductance  between  node  (i,j,k)  and  node  (i,j,k+l)  is  read  as 
CV(i,j,k).  There  are  three  other  conductances  connected  to  node  (i,j,k),  but  these  are  associated 
with  neighboring  nodes.  For  example,  the  conductance  between  node  (i,j,k)  and  node  (i,j-l,k)  is 
CC(i,j-l,k).  CDTR  and  CDTC  values  are  indexed  the  same  as  CR  and  CC  respectively. 
CDTR(i,j,k)  is  CDT  between  nodes  (i,j,k)  and  (i,j+l,k);  CDTC(i,j,k)  is  CDT  between  nodes  (i,j,k) 
and  (i+l,j,k). 

Although  conductance  and  CDT  values  are  indexed  according  to  node  location,  it  is  important  to 
remember  that  the  values  represent  the  material  between  nodes.  There  is  no  need  for  vertical 
conductance  values  for  the  bottom  layer  because  the  vertical  conductance  assigned  to  the  nodes  in 
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the  bottom  layer  represents  conductance  between  the  nodes  in  that  layer  and  a  nonexistent  layer 
below.  Accordingly,  the  vertical  conductance  for  the  bottom  layer  is  not  read;  CV  for  the  bottom 
layer  must  not  be  included  in  the  input  data.  For  each  layer,  either  CR  or  CDTR  must  be  read; 
however,  the  values  for  the  last  column  are  not  used.  That  is,  values  for  nodes  (i,NCOL,k),  in 
which  NCOL  is  the  maximum  column  number,  are  not  used.  This  is  because  these  values  represent 
conductance  between  nodes  in  the  last  column  and  the  next  column,  which  does  not  exist  because 
it  is  outside  of  the  model.  Although  values  must  be  specified  for  the  last  column,  the  values  for  the 
last  column  are  ignored.  Similarly  for  CC  and  CDTC,  the  values  for  the  last  row,  that  is  for  nodes 
(NROW,j,k)  in  which  NROW  is  the  maximum  row  number,  must  be  included  as  part  of  the  input 
array;  however,  these  values  are  ignored  because  they  represent  conductance  to  nodes  that  are 
outside  of  the  model.  A  no-flow  boundary  is  automatically  imposed  around  the  entire  model 
perimeter. 

Storage  capacity  is  not  needed  for  steady-state  simulations.  A  flag  (ISS)  has  been  provided  to 
allow  the  user  to  specify  that  a  simulation  is  steady  state.  When  ISS  is  nonzero,  space  for  storage 
capacity  is  not  allocated,  storage  capacity  is  not  read,  and  storage  calculations  are  skipped. 

Four  values  for  the  layer-type  code  are  accepted  by  the  GFD  Package.  These  values  identify  which 
of  the  various  combinations  of  the  options  provided  by  the  package  are  to  be  used  for  a  layer.  The 
layer- type  codes  are  stored  in  the  one-dimensional  array  LAYCON.  The  code  values  and  the 
corresponding  options  invoked  are  given  below.  The  layer-type  code  is  identical  in  meaning  to  the 
layer-type  code  in  the  BCF  Package. 

Layer-type  code  0  indicates  that  horizontal  conductance  and  storage  capacity  are  constant.  There 
is  no  provision  for  adjusting  conductance,  limiting  vertical  flow  from  above,  or  adjusting  storage 
terms  as  water  level  varies.  This  layer  type  is  usually  used  to  simulate  confined  conditions.  This 
option  could  also  be  used  to  simulate  unconfined  conditions  provided  that  simulated  water-level 
changes  are  a  small  fraction  of  saturated  thickness  and  flow  from  an  overlying  layer  (if  present) 
can  be  assumed  to  be  negligible. 

Input  required  for  each  node  in  the  layer: 

SCI  -  Confined  or  unconfined  storage  capacity  (only  for  transient 
simulations), 

CR  —  Conductance  between  each  node  and  the  next  adjacent  node  along  a 
row, 

CC  —  Conductance  between  each  node  and  the  next  adjacent  node  along  a 
column,  and 

CV  —  Conductance  between  each  node  in  the  layer  and  each  node  in  the  layer 
below  (only  if  there  is  a  layer  below). 

Layer-type  code  1  indicates  a  strictly  unconfined  layer,  and  can  be  used  only  for  the  top  layer  of  a 
model.  Horizontal  conductance  is  calculated  each  iteration  as  the  product  of  saturated  thickness 
and  conductance  divided  by  saturated  thickness;  however,  there  is  no  upper  limit  on  saturated 
thickness,  and  there  is  no  provision  for  storage  term  conversion.  Also,  there  is  never  any  limitation 
on  vertical  flow  from  layers  above  because  this  layer  type  is  valid  only  for  the  top  layer. 
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Input  required  for  each  node  in  the  layer: 

SCI  --  Unconfined  storage  capacity  (only  for  transient  simulations), 

CDTR  —  Conductance  divided  by  thickness  between  each  node  and  the  next 
adjacent  node  along  a  row, 

CDTC  --  Conductance  divided  by  thickness  between  each  node  and  the  next 
adjacent  node  along  a  column, 

BOT  —  Elevation  of  the  bottom  of  the  aquifer,  and 

CV  —  Conductance  between  each  node  in  the  layer  and  each  node  in  the  layer 
below  (only  if  there  is  a  layer  below). 

Layer-type  code  2  indicates  a  layer  that  may  alternate  between  confined  and  unconfined,  but 
conductance  is  held  constant  throughout  the  simulation.  This  is  useful  in  situations  where  storage 
term  conversion  and  limitation  of  flow  from  above  under  dewatered  conditions  are  desirable,  and 
yet  saturated  thickness  changes  are  small  compared  to  the  saturated  thickness  so  that  recalculation 
of  conductance  as  a  function  of  saturated  thickness  changes  is  not  necessary. 

Input  required  for  each  node  in  the  layer. 

SCI  —  Storage  capacity  to  be  used  when  nodes  are  confined  (only  for  transient 
simulations), 

CR  —  Conductance  between  each  node  and  the  next  adjacent  node  along  a 
row, 

CC  —  Conductance  between  each  node  and  the  next  adjacent  node  along  a 
column, 

SC2  -  Storage  capacity  to  be  used  when  nodes  are  unconfined  (only  for 
transient  simulations), 

TOP  -  Elevation  of  the  top  of  the  aquifer,  and 

CV  --  Conductance  between  each  node  in  the  layer  and  each  node  in  the  layer 
below  (only  if  there  is  a  layer  below). 

Layer-type  code  3  indicates  a  layer  that  may  alternate  between  confined  and  unconfined,  and  all  of 
the  options  associated  with  water-table  conditions  are  incorporated.  Conductance  is  calculated 
each  iteration  as  the  product  of  saturated  thickness  and  conductance  divided  by  saturated  thickness. 
Storage  terms  can  convert  back  and  forth  between  confined  and  unconfined,  and  vertical  flow  from 
above  is  limited  at  nodes  that  are  unconfined. 

Input  required  for  each  node  in  the  layer: 

SCI  —  Storage  capacity  to  be  used  when  nodes  are  confined  (only  for  transient 
simulations), 

CDTR  —  Conductance  divided  by  thickness  between  each  node  and  the  next 
adjacent  node  along  a  row, 

CDTC  -  Conductance  divided  by  thickness  between  each  node  and  the  next 
adjacent  node  along  a  column, 

BOT  -  Elevation  of  the  bottom  of  the  aquifer, 

SC2  —  Storage  capacity  to  be  used  when  nodes  are  unconfined  (only  for 
transient  simulations), 

TOP  —  Elevation  of  the  top  of  the  aquifer,  and 
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CV  -  Conductance  between  each  node  in  the  layer  and  each  node  in  the  layer 
below  (only  if  there  is  a  layer  below). 

Only  the  required  arrays  as  indicated  by  the  LAYCON  value  for  each  layer  are  read.  Care  must  be 
taken  to  ensure  that  only  the  required  data  is  included  in  the  input  data.  Otherwise,  the  data  will 
be  incorrectly  read  by  the  program. 

The  four  layer-type  designations  allow  the  user  to  control  the  complexity  of  the  formulations. 
Layer-type  codes  1-3  result  in  greater  usage  of  computer  memory  and  execution  time  than  does 
layer-type  code  0.  Layer-type  code  3  uses  the  most  memory  and  execution  time.  Also,  layer-type 
codes  1-3  can  cause  numerical  instability  that  can  prevent  the  iterative  solvers  from  converging  to 
a  solution,  and  therefore,  layer-type  codes  1-3  should  be  used  only  when  necessary. 


INPUT  INSTRUCTIONS 

The  input  to  the  General  Finite-Difference  Flow  (GFD)  Package  is  listed  by  “items”  of  input.  An 
item  defines  one  or  more  input  parameters,  and  each  item  requires  one  or  more  input  records.  An 
input  parameter  may  be  an  individual  value,  or  it  may  be  an  array  having  multiple  values.  Each 
item  is  numbered.  The  description  of  an  item  contains  a  line  that  names  the  parameters  and  a  line 
showing  the  format  of  the  parameters.  For  parameters  that  are  individual  values  and  some  of  the 
parameters  that  are  arrays,  the  format  is  given  in  standard  FORTRAN  notation.  The  format  for  the 
remainder  of  the  arrays  is  defined  by  naming  the  utility  module  that  reads  them.  Details  about  the 
array  reading  utility  modules  are  contained  in  McDonald  andHarbaugh  (1988,  p.  14-4).  Following 
the  list  of  items  is  a  list  of  definitions  of  input  parameters.  The  variables  and  arrays  used  in  the 
program  have  the  same  names  as  the  parameters  they  represent. 

Input  for  the  GFD  Package  is  read  from  the  unit  specified  in  IUNIT(14).  IUNIT  is  an  array  of  file 
unit  numbers  that  is  read  by  the  Basic  Package.  A  positive  number  for  IUNIT(14)  indicates  that 
GFD  is  to  be  used,  and  that  number  is  the  unit  number  on  which  the  input  data  are  to  be  read.  See 
McDonald  and  Harbaugh  (1988,  p.  3-25)  for  an  explanation  of  the  IUNIT  array.  Note  that  a 
different  element  within  IUNIT  may  be  used  if  IUNIT(14)  is  already  being  used  by  another 
package.  (See  Module  Documentation  section.)  When  executing  the  model,  either  the  GFD  or  the 
BCF  Package  must  be  used.  IUNIT(l)  is  the  input  unit  number  for  the  BCF  Package.  Thus,  either 
IUNIT(l)  or  IUNIT(14)  must  be  positive,  but  not  both.  The  data  are  read  once  at  the  start  of  a 
simulation. 
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Module  GFD1 AL  reads  items  1-2: 


1.  Data:  ISS  IGFDCB 

Format:  110  110 

2.  Data:  LAYCON(NLAY)  (Maximum  of  80  layers) 

Format:  4012 

(If  there  are  40  or  fewer  layers,  use  one  record;  otherwise,  use  two  records.) 
Module  GFD1RP  reads  items  3-13: 

3.  Data:  DELR(NCOL) 

Module:  U1DREL 

4.  Data:  DELC(NROW) 

Module:  U1DREL 

A  subset  of  the  following  two-dimensional  arrays  are  used  to  describe  each  layer.  The  arrays 
needed  for  each  layer  depend  on  the  layer  type  code  (LAYCON)  and  whether  the  simulation  is 
transient  (ISS  =  0)  or  steady  state  (ISS  not  0).  If  an  array  is  not  needed,  it  must  be  omitted.  All  of 
the  arrays  (items  5-13)  for  layer  1  are  read  first;  then  all  of  the  arrays  for  layer  2,  etc. 

IF  THE  SIMULATION  IS  TRANSIENT  (ISS  not  0) 

5.  Data:  SCl(NCOL,NROW) 

Module:  U2DREL 

IF  THE  LAYER-TYPE  CODE  (LAYCON)  IS  0  OR  2 

6.  Data:  CR(NCOL,NROW) 

Module:  U2DREL 

7.  Data:  CC(NCOL,NROW) 

Module:  U2DREL 

IF  THE  LAYER-TYPE  CODE  (LAYCON)  IS  1  OR  3 

8.  Data:  CDTR(NCOL,NROW) 

Module:  U2DREL 

9.  Data:  CDTC(NCOL,NROW) 

Module:  U2DREL 

10.  Data:  BOT(NCOL,NROW) 

Module:  U2DREL 
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IF  THE  LAYER  IS  NOT  THE  BOTTOM  LAYER 


11.  Data:  CV(NCOL,NROW) 

Module:  U2DREL 

IF  THE  SIMULATION  IS  TRANSIENT  (ISS  not  0)  AND  THE  LAYER-TYPE  CODE  (LAYCON) 
IS  2  OR  3 

12.  Data:  SC2(NCOL,NROW) 

Module:  U2DREL 

IF  THE  LAYER-TYPE  CODE  IS  2  OR  3 

13.  Data:  TOP(NCOL,NROW) 

Module:  U2DREL 

Definition  of  Input  Parameters 

ISS — is  the  steady-state  flag. 

If  ISS  is  not  0,  the  simulation  is  steady  state. 

If  ISS  =  0,  the  simulation  is  transient. 

IGFDCB-is  a  flag  and  a  unit  number. 

If  IGFDCB  >  0,  it  is  the  unit  number  on  which  cell-by-cell  flow  terms  will  be 

recorded  whenever  ICBCFL  is  set.  ICBCFL  is  an  input  item  in  the 
Output  Control  option  of  the  model  (McDonald  and  Harbaugh, 
1988,  p.4-14). 

If  IGFDCB  =  0,  cell-by-cell  flow  terms  will  not  be  printed  or  recorded. 

If  IGFDCB  <  0,  flow  for  each  constant- head  cell  will  be  printed  whenever  ICBCFL 
is  set. 

LAYCON— is  the  layer- type  code  array.  Each  element  holds  the  code  for  the  respective  layer. 

Read  one  value  for  each  layer.  There  is  a  limit  of  80  layers.  Leave  unused  elements 
blank. 

0  -  CONFINED -Horizontal  conductance  and  storage  capacity  of  the  layer  are  constant 
for  the  entire  simulation. 

1  -  UNCONFINED-Horizontal  conductance  of  the  layer  varies.  It  is  calculated  from 
the  saturated  thickness  and  the  conductance  divided  by 
thickness.  The  storage  capacity  is  constant;  this  value  is  valid 
only  for  layer  1 . 
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2  -  CONFINED/UNCONFINED— Horizontal  conductance  of  the  layer  is  constant.  The 

storage  capacity  may  be  either  confined  or 
unconfined.  Vertical  flow  from  above  is  limited 
while  unconfined. 

3  -  CONFINED/UNCONFINED-Horizontal  conductance  of  the  layer  varies.  It  is 

calculated  from  the  saturated  thickness  and 
conductance  divided  by  thickness.  The  storage 
capacity  may  be  either  confined  or  unconfined. 
Vertical  flow  from  above  is  limited  while 
unconfined. 

DELR-is  the  cell  width  along  rows.  Read  one  value  for  each  of  the  NCOL  columns. 

DELC-is  the  cell  width  along  columns.  Read  one  value  for  each  of  the  NROW  rows. 

SCI— is  the  primary  storage  capacity.  Read  SCI  only  for  a  transient  simulation  (steady-state  flag, 
ISS,  is  0).  If  LAYCON  is  1  (strictly  unconfined),  SCI  is  equal  to  specific  yield  multiplied 
by  cell  area;  if  LAYCON  is  2  or  3,  SCI  is  confined  storage  coefficient  multiplied  by  cell  area. 
If  LAYCON  is  0,  the  layer  is  normally  confined  and  SCI  would  then  be  the  confined  storage 
coefficient  multiplied  by  cell  area.  However,  layer-type  0  can  be  used  to  simulate  an 
unconfined  aquifer  if  saturated  thickness  changes  are  small  compared  to  total  saturated 
thickness  and  if  there  is  negligible  flow  from  above;  in  this  case,  SCI  would  be  specific  yield 
multiplied  by  cell  area. 

CR-is  the  conductance  along  rows.  Read  CR  only  for  layers  where  LAYCON  is  0  or  2.  The  value 
of  CR  for  node  (i,j,k)  represents  conductance  between  that  node  and  node  (i,j+l,k). 

CC-is  the  conductance  along  columns.  Read  CC  only  for  layers  where  LAYCON  is  0  or  2.  The 
value  of  CC  for  node  (i,j,k)  represents  conductance  between  that  node  and  node  (i+l,j,k). 

CDTR— is  the  conductance  divided  by  thickness  along  rows.  Read  CDTR  only  for  layers  where 
LAYCON  is  1  or  3.  The  value  of  CDTR  for  node  (i,j,k)  represents  conductance  divided  by 
thickness  between  that  node  and  node  (i,j+l,k). 

CDTC— is  the  conductance  divided  by  thickness  along  columns.  Read  CDTC  only  for  layers  where 
LAYCON  is  1  or  3.  The  value  of  CDTC  for  node  (i,j,k)  represents  conductance  divided  by 
thickness  between  that  node  and  node  (i+l,j,k). 

BOT— is  the  elevation  of  the  aquifer  bottom.  Read  BOT  only  for  layers  where  LAYCON  is  1  or  3. 

CV— is  the  vertical  conductance.  The  value  of  CV  for  node  (i,j,k)  represents  conductance  between 
that  node  and  node  (i,j,k+l).  Because  there  is  not  a  layer  beneath  the  bottom  layer,  CV  is  not 
required  for  the  bottom  layer  and  must  not  be  included  in  the  data. 

SC2— is  the  secondary  storage  capacity.  Read  SC2  only  for  layers  where  LAYCON  is  2  or  3  and 
only  if  a  transient  simulation  (steady-state  flag,  ISS,  is  0).  The  secondary  storage  capacity  is 
always  specific  yield  multiplied  by  cell  area. 

TOP-is  the  elevation  of  the  aquifer  top.  Read  TOP  only  for  layers  where  LAYCON  is  2  or  3. 
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SAMPLE  PROBLEM 


To  demonstrate  the  use  of  the  GFD  Package,  the  sample  problem  from  McDonald  and  Harbaugh 
(1988,  Appendix  D)  has  been  executed  using  the  GFD  Package.  The  GFD  input  data  is  shown  in 
Table  1.  Layer  1,  simulated  as  an  unconfined  aquifer  (LAYCON=l),  requires  CDT  values,  which 
are  calculated  as  hydraulic  conductivity  times  cell  width  divided  by  distance  between  nodes 
(equation  3).  Horizontal  cell  dimensions  are  a  constant  5000  ft  in  both  horizontal  directions,  and 
assuming  that  nodes  are  still  centered  in  cells,  the  distance  between  adjacent  nodes  is  5000  ft. 
Thus,  cell  width  divided  by  distance  between  nodes  is  one  because  of  the  uniform  cell  spacing  in 
both  directions.  The  resulting  CDT  values  are  uniform  in  both  directions  for  all  cells  in  the  layer; 
CDTR  and  CDTC  for  layer  1  are  .001  ft/s.  Layers  2  and  3  are  simulated  as  confined  aquifers,  and 
conductance  is  required.  For  each  of  these  layers,  transmissivity  is  uniform.  Conductance  is 
calculated  as  the  product  of  transmissivity  times  cell  width  divided  by  distance  between  nodes 
(equation  2).  The  result  is  that  conductances  along  rows  (CR)  and  columns  (CC)  are  the  same  for 
all  cells  in  each  layer.  Conductance  in  layer  2  is  0.01  ft?/s,  and  0.02  ft2/s  for  layer  3.  Vertical 
conductance  (CV)  is  the  value  for  VCONT  from  the  BCF  Package  multiplied  by  cell  area.  The 
resulting  value  for  CV  between  layers  1  and  2  is  0.5  ft2/s;  the  value  is  0.25  ft2/s  between  layers  2 
and  3. 


Table  1.  —  Sample  problem  input  data  for  the  GFD  Package 

[Note  that  parameter  names,  such  as  DELC,  are  included  in  the  following 
records  in  order  to  help  identify  the  data.  Such  labels,  or  any  desired 
information,  can  be  included  in  the  actual  input  data  read  by  the  model 
program  provided  that  the  added  information  does  not  occupy  any  area  of  a 
record  used  for  data.  Extra  records  cannot  be  included.] 


1 

10  0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

5000  . 
5000  . 
.001 
.001 
-150  . 
.5 
.01 
.01 
.25 
.02 


0  .02 


ISS, IGFDBD 
LAYCON 


DELR 

DELC 

CDTR,  layer  1 
CDTC,  layer  1 
BOT,  layer  1 
CV,  layer  1 
CR,  layer  2 
CC,  layer  2 
CV,  layer  2 
CR,  layer  3 
CC,  layer  3 


The  input  to  the  Basic  Package  is  unchanged  except  IUNIT(l)  is  set  to  0  and  IUNIT(14)  is  set  to 
1 1  to  deactivate  the  BCF  Package  and  activate  the  GFD  Package.  Input  for  all  the  other  packages 
is  unchanged  from  that  shown  in  McDonald  and  Harbaugh  (1988,  Appendix  D). 


The  results  from  running  the  sample  problem  using  the  GFD  Package  are  shown  in  Table  2.  The 
heads  are  nearly  equal  to  those  from  McDonald  and  Harbaugh  (1988);  the  differences  can  be 
attributed  to  the  different  formulation  of  conductance  for  the  top  layer.  The  GFD  Package  assumes 
a  uniform  CDT  and  linearly  varying  saturated  thickness  between  nodes.  The  BCF  Package 
assumes  a  constant  saturated  thickness  and  hydraulic  conductivity  in  each  cell.  If  heads  were 
exacdy  the  same  in  both  simulations,  the  assumptions  used  in  the  GFD  Package  would  result  in 
slightly  higher  conductances  than  the  conductances  calculated  by  the  BCF  Package.  The  higher 
conductances  would  cause  too  much  water  to  flow  through  the  system,  which  would  cause  heads 
to  decrease  in  order  to  correctly  solve  the  flow  equation.  Thus,  heads  in  the  GFD  simulation  are 
slightly  lower  than  heads  in  the  BCF  simulation. 
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Table  2. —Output  from  the  sample  problem 


U.S.  GEOLOGICAL  SURVEY  MODULAR  FINITE-DIFFERENCE  GROUND-WATER  MODEL 

SAMPLE  PROBLEM  FROM  TWRI  6-Al  —  HEAD  DIFFERS  SLIGHTLY  DUE  TO  MODIFIED  CONDUCTANCE  CALCULATION  FOR  WATER-TABLE  LAYER 
3  LAYERS  15  ROWS  15  COLUMNS 

1  STRESS  PERIOD (S)  IN  SIMULATION 
MODEL  TIME  UNIT  IS  SECONDS 

I/O  UNITS: 

ELEMENT  OF  IUNIT:  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24 
I/O  UNIT:  0  12  13  0000  18  19  00000  11  000000000 

BAS1  —  BASIC  MODEL  PACKAGE,  VERSION  1,  9/1/87  INPUT  READ  FROM  UNIT  1 
ARRAYS  RHS  AND  BUFF  WILL  SHARE  MEMORY. 

START  HEAD  WILL  NOT  BE  SAVED  —  DRAWDOWN  CANNOT  BE  CALCULATED 
5892  ELEMENTS  IN  X  ARRAY  ARE  USED  BY  BAS 
5892  ELEMENTS  OF  X  ARRAY  USED  OUT  OF  350000 

GFD1  —  GENERAL  FINITE-DIFFERENCE  FLOW  PACKAGE,  VERSION  1,  9/19/89  INPUT  READ  FROM  UNIT  11 
STEADY-STATE  SIMULATION 
LAYER  AQUIFER  TYPE 


1  1 

2  0 

3  0 

675  ELEMENTS  IN  X  ARRAY  ARE  USED  BY  GFD 
6567  ELEMENTS  OF  X  ARRAY  USED  OUT  OF  350000 

WEL1  —  WELL  PACKAGE,  VERSION  1,  9/1/87  INPUT  READ  FROM  12 

MAXIMUM  OF  15  WELLS 

60  ELEMENTS  IN  X  ARRAY  ARE  USED  FOR  WELLS 

6627  ELEMENTS  OF  X  ARRAY  USED  OUT  OF  350000 

DRNl  —  DRAIN  PACKAGE,  VERSION  1,  9/1/87  INPUT  READ  FROM  UNIT  13 

MAXIMUM  OF  9  DRAINS 

45  ELEMENTS  IN  X  ARRAY  ARE  USED  FOR  DRAINS 

6672  ELEMENTS  OF  X  ARRAY  USED  OUT  OF  350000 

RCH1  —  RECHARGE  PACKAGE,  VERSION  1,  9/1/87  INPUT  READ  FROM  UNIT  18 
OPTION  1  —  RECHARGE  TO  TOP  LAYER 

225  ELEMENTS  OF  X  ARRAY  USED  FOR  RECHARGE 

6897  ELEMENTS  OF  X  ARRAY  USED  OUT  OF  350000 

SIP1  —  STRONGLY  IMPLICIT  PROCEDURE  SOLUTION  PACKAGE,  VERSION  1,  9/1/87  INPUT  READ  FROM  UNIT  19 
MAXIMUM  OF  50  ITERATIONS  ALLOWED  FOR  CLOSURE 
5  ITERATION  PARAMETERS 

2905  ELEMENTS  IN  X  ARRAY  ARE  USED  BY  SIP 

9802  ELEMENTS  OF  X  ARRAY  USED  OUT  OF  350000 


Table  2.~0utput  from  the  sample  problem  (Cont.) 


SAMPLE  PROBLEM  FROM  TWRI  6-A1  —  HEAD  DIFFERS  SLIGHTLY  DUE  TO  MODIFIED  CONDUCTANCE  CALCULATION  FOR  WATER-TABLE  LAYER 

BOUNDARY  ARRAY  FOR  LAYER  1  WILL  BE  READ  ON  UNIT  1  USING  FORMAT:  (1513) 

1  2  3  4  5  6  7  8  9  10  11  12  13  14  15 

1-111111111111111 

2-111111111111111 

3- 111111111111111 

4- 111111111111111 

5- 111111111111111 

6- 111111111111111 

7- 111111111111111 

8- 111111111111111 
9-111111111111111 

10  -111111111111111 

11  -111111111111111 
12  -11  1111111111111 

13  -111111111111111 

14  -111111111111111 

15  -111111111111111 


Table  2.~< 


1  2  3  4  5  6 

1-111111 

2-111111 

3- 111111 

4- 111111 

5- 111111 

6- 111111 

7- 111111 

8- 111111 

9-111111 

10  -1  1  1  1  1  1 

11  -1  1  1  1  1  1 

12  -1  1  1  1  1  1 

13  -1  1  1  1  1  1 

14  -1  1  1  1  1  1 

15  -1  1  1  1  1  1 


BOUNDARY  ARRAY  FOR  LA' 

7  8  9  10  11  12  13  14 

11111111 

11111111 

11111111 

11111111 

11111111 

11111111 

11111111 

11111111 

11111111 

11111111 

11111111 

11111111 

11111111 

11111111 

11111111 


15 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 


Table  2.-0utput  from  the  sample  problem  (Cont.) 


AQUIFER  HEAD  WILL  BE  SET  TO  999.99 


BOUNDARY  ARRAY  = 
AT  ALL  NO-FLOW  NODES  (IBOUND=0) . 


1  FOR  LAYER  3 


INITIAL 

HEAD  = 

0.0000000 

FOR 

LAYER 

1 

INITIAL 

HEAD  = 

0.0000000 

FOR 

LAYER 

2 

INITIAL 

HEAD  = 

0.0000000 

FOR 

LAYER 

3 

DEFAULT  OUTPUT  CONTROL  —  THE  FOLLOWING  OUTPUT  COMES  AT  THE  END  OF  EACH  STRESS  PERIOD: 
TOTAL  VOLUMETRIC  BUDGET 
HEAD 


DELR 

= 

5000.000 

DELC 

= 

5000.000 

COND/THICK 

ALONG 

ROWS 

= 

0 . 9999999E-03 

FOR 

LAYER 

1 

COND/THICK 

ALONG 

COLS 

= 

0 . 9999999E-03 

FOR 

LAYER 

1 

BOTTOM 

= 

-150.0000 

FOR 

LAYER 

1 

VERTICAL  * 

CONDUC' 

TANCF. 

= 

0.5000000 

FOR 

LAYER 

1 

CONDUCTANCE 

ALONG 

ROWS 

= 

0 .1000000E-01 

FOR 

LAYER 

2 

CONDUCTANCE 

ALONG 

COLS 

= 

0 . 1000000E-01 

FOR 

LAYER 

2 

VERTICAL  ' 

CONDUCTANCE 

= 

0.2500000 

FOR 

LAYER 

2 

CONDUCTANCE 

ALONG 

ROWS 

= 

0 .2000000E-01 

FOR 

LAYER 

3 

CONDUCTANCE 

ALONG 

COLS 

= 

0 .2000000E-01 

FOR 

LAYER 

3 

SOLUTION  BY  THE  STRONGLY  IMPLICIT  PROCEDURE 


MAXIMUM  ITERATIONS  ALLOWED  FOR  CLOSURE  =  50 

ACCELERATION  PARAMETER  =  1.0000 

HEAD  CHANGE  CRITERION  FOR  CLOSURE  =  0.10000E-02 

SIP  HEAD  CHANGE  PRINTOUT  INTERVAL  =  1 


Table  2.— Output  from  the  sample  problem  (Cont.) 


5  ITERATION  PARAMETERS  CALCULATED  FROM  SPECIFIED  WSEED  =  0.00100000  : 

0.0000000E  +  00  0 . 8221720E  +  00  0 . 9683772E  +  00  0 . 9943765E+00  0 . 9990000E  +  00 


STRESS  PERIOD  NO.  1,  LENGTH  =  86400.00 


NUMBER  OF  TIME  STEPS  =  1 

MULTIPLIER  FOR  DELT  =  1.000 

INITIAL  TIME  STEP  SIZE  =  86400.00 


LAYER 

ROW 

COL 

STRESS  RATE 

WELL  ] 

3 

5 

11 

-5.0000 

1 

2 

4 

6 

-5.0000 

2 

2 

6 

12 

-5.0000 

3 

1 

9 

8 

-5.0000 

4 

1 

9 

10 

-5.0000 

5 

1 

9 

12 

-5.0000 

6 

1 

9 

14 

-5.0000 

7 

1 

11 

8 

-5 . 0000 

8 

1 

11 

10 

-5.0000 

9 

1 

11 

12 

-5.0000 

10 

1 

11 

14 

-5 . 0000 

11 

1 

13 

8 

-5.0000 

12 

1 

13 

10 

-5.0000 

13 

1 

13 

12 

-5.0000 

14 

1 

13 

14 

-5.0000 

15 

9  DRAINS 


LAYER 

ROW 

COL 

ELEVATION 

CONDUCTANCE 

DRAIN  NO. 

1 

8 

2 

0.0000 

1.000 

1 

1 

8 

3 

0.0000 

1.000 

2 

1 

8 

4 

10.00 

1.000 

3 

1 

8 

5 

20.00 

1.000 

4 

1 

8 

6 

30.00 

1.000 

5 

1 

8 

7 

50.00 

1.000 

6 

1 

8 

8 

70.00 

1.000 

7 

1 

8 

9 

90.00 

1.000 

8 

1 

8 

10 

100.0 

1.000 

9 

RECHARGE 


0 . 3000000E-07 


Table  2.— Output  from  the  sample  problem  (Cont.) 


31  ITERATIONS  FOR  TIME  STEP  3  IN  STRESS  PERIOD  1 


MAXIMUM  HEAD  CHANGE  FOR  EACH  ITERATION: 


HEAD  CHANGE 

LAYER, 

.  ROW, 

COL 

HEAD  CHANGE 

LAYER, 

.ROW, 

COL 

HEAD  CHANGE 

LAYER, 

.  ROW, 

.COL 

HEAD  CHANGE 

LAYER, 

,  ROW, 

COL 

HEAD  CHANGE 

LAYER, 

.ROW, 

COL 

-22.41 

(  3, 

5, 

11) 

12.48 

(  1, 

1, 

15) 

13.39 

(  3, 

1, 

14) 

48.20 

(  1, 

1, 

15) 

35.87 

(  3, 

1, 

13) 

2.481 

(  1, 

9, 

14) 

1.429 

(  3, 

10, 

13) 

6.210 

(  1, 

12, 

14) 

7.405 

(  3, 

11, 

14) 

13.64 

(  1, 

15, 

15) 

0.5474 

(  3, 

8, 

7) 

0.4805 

(  2, 

6, 

9) 

0.4700 

<  3, 

5, 

10) 

2.015 

(  1, 

11, 

14) 

2.296 

(  3, 

5, 

13) 

0.1105 

(  1, 

13, 

12) 

0 . 7037E-01 

(  3, 

12, 

11) 

0.2809 

(  1, 

14, 

14) 

0.3130 

(  3, 

13, 

14) 

0.3299 

(  1, 

15, 

15) 

0 . 7  80  9E-02 

<  1, 

13, 

12) 

0 . 1580E-01 

(  2, 

11, 

ID 

0 . 1770E-01 

{  3, 

11, 

10) 

0 .7875E-01 

(  1, 

14, 

14) 

0 . 8456E-01 

(  3, 

7, 

14) 

0.4147E-02 

(  1, 

13, 

14) 

0 .2542E-02 

(  3, 

14, 

15) 

0 . 9710E-02 

<  1, 

14, 

14) 

0.107  6E-01 

(  3, 

13, 

14) 

0.1019E-01 

(  1, 

15, 

15) 

0 . 2408E-03  (  1,  13,  12) 


Table  2.— Output  from  the  sample  problem  (Cont.) 


HEAD 

IN  LAYER 

1  AT  END  OF  TIME 

STEP  1 

IN  STRESS  PERIOD 

1 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

1 

0.0000 

24.86 

43.91 

59.16 

71.74 

82.44 

91.84 

99.97 

106.9 

112.6 

117.3 

121.2 

124.2 

126.3 

127.4 

2 

0.0000 

24.36 

43.01 

57.89 

70.09 

80.50 

90.05 

98.33 

105.2 

110.9 

115.6 

119.5 

122.7 

124 . 9 

126.0 

3 

0.0000 

23.37 

41.21 

55.35 

66.70 

76.14 

86.45 

95.14 

102.1 

107.5 

112.0 

116.1 

119.6 

122.0 

123.3 

4 

0.0000 

21.85 

38.53 

51.68 

61.73 

67.97 

81.28 

90.69 

97.58 

102.5 

106.0 

110.7 

114.9 

117.8 

119.3 

5 

0.0000 

19.68 

34.85 

47.25 

57.63 

66.68 

77.03 

85.71 

92.17 

96.09 

97.24 

103.1 

108.8 

112.5 

114.3 

6 

0.0000 

16.47 

29.45 

40.85 

51.25 

61.16 

71.14 

79.80 

86.42 

90.77 

92.98 

94.18 

102.0 

106.4 

108.4 

7 

0.0000 

11.53 

21.07 

31.17 

41.36 

51.81 

63.04 

72.64 

79.91 

84.87 

88.55 

91.62 

96.38 

99.77 

101.8 

8 

0.0000 

3.479 

6.830 

16.25 

26.30 

36.96 

52.57 

64.27 

72.48 

77.21 

81.94 

84.96 

89.22 

91.67 

94.29 

9 

0.0000 

10.53 

19.08 

28.10 

36.90 

45.25 

52.93 

55.35 

65.11 

66.03 

73.89 

73.74 

80.79 

80.12 

86.45 

10 

0.0000 

14.60 

25.83 

35.35 

43.46 

50.08 

54 . 90 

57.51 

62.91 

65.51 

70.34 

72.39 

76.67 

78.21 

81.75 

11 

0.0000 

17.08 

29.91 

39.97 

47.74 

53.20 

55.77 

53.29 

60.23 

59.25 

66.39 

65.41 

72.18 

71.00 

77.58 

12 

0.0000 

18.64 

32.51 

43.02 

50.77 

55.88 

58.29 

58.43 

61.89 

63.14 

67.08 

68.46 

72.25 

73.41 

76.81 

13 

0.0000 

19.62 

34.18 

45.08 

52 . 96 

57.99 

59.86 

56.71 

62.55 

60.87 

67.18 

65.71 

71.86 

70.31 

76.44 

14 

0.0000 

20.22 

35.21 

46.42 

54.56 

60.03 

63.13 

64.47 

67.21 

68.74 

71.60 

73.13 

75.79 

76.98 

79.05 

15 

0.0000 

20.50 

35.71 

47.10 

55.43 

61.21 

64 . 98 

67.48 

69.90 

71.97 

74.24 

76.17 

78.17 

79.62 

80.78 

Table  2. -Output  from  the  sample  problem  (Cont.) 


HEAD 

IN  LAYER  2 

AT  END  OF  TIME 

STEP  1 

IN  STRESS  PERIOD 

1 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

1 

0.0000 

24.58 

43.63 

58.92 

71.53 

82.24 

91.65 

99.80 

106.7 

112.4 

117.2 

121.1 

124.1 

126.2 

127.2 

2 

0.0000 

24.09 

42.73 

57.65 

69.87 

80.29 

89.86 

98.16 

105.1 

110.8 

115.4 

119.4 

122.5 

124.7 

125.9 

3 

0.0000 

23 . 10 

40.94 

55.11 

66.46 

75.70 

86.22 

94.96 

101 . 9 

107.4 

111.8 

115.9 

119.4 

121 . 9 

123.1 

4 

0.0000 

21.59 

38.26 

51.43 

61.28 

60.11 

80.84 

90.49 

97.39 

102.2 

105.3 

110.4 

114 . 7 

117.7 

119.2 

5 

0.0000 

19.43 

34.59 

47.01 

57.38 

66.24 

76.80 

85.51 

91.94 

95.36 

91.04 

102.1 

108.5 

112.3 

114.1 

6 

0.0000 

16.23 

29.20 

40.60 

51.02 

60.93 

70.93 

79.60 

86.23 

90.49 

92.01 

86.18 

101.6 

106.2 

108.2 

7 

0.0000 

11.36 

20.92 

31.01 

41.21 

51.66 

62.86 

72.44 

79.72 

CO 

00 

88.30 

91.19 

96.17 

99.60 

101.6 

8 

0.0000 

4.205 

8.326 

17.58 

27.57 

38.24 

52.93 

64.15 

72.30 

77.07 

81.76 

84.81 

89.05 

91.54 

94.12 

9 

0.0000 

10.37 

18 . 94 

27.96 

36.77 

45.14 

52.84 

56.10 

65.04 

66.75 

73.82 

74.44 

80.72 

80.80 

86.33 

10 

0.0000 

14.37 

25.58 

35.12 

43.24 

49.88 

54.73 

57.45 

62.76 

65.44 

70.20 

72.33 

76.53 

78.15 

81.59 

11 

0.0000 

16.84 

29.66 

39.74 

47.52 

53.01 

55.64 

54.06 

60.16 

60.00 

66.33 

66.14 

72.12 

71.71 

77.47 

12 

0.0000 

18.39 

32.26 

42.80 

50.56 

55.69 

58.12 

58.37 

61.74 

63.08 

66.93 

68.40 

72.11 

73.36 

76.65 

13 

0.0000 

19.38 

33 . 93 

44.86 

52.75 

57.80 

59.74 

57.46 

62.48 

61.61 

67.12 

66.44 

71.80 

71.01 

76.33 

14 

0.0000 

19.97 

34.96 

46.20 

54.35 

59.83 

62.95 

64.35 

67.04 

68.62 

71.43 

73.01 

75.63 

76.87 

78.88 

15 

0.0000 

20.25 

35.46 

46.88 

55.22 

61.02 

64.79 

67.30 

69.72 

71.80 

74.07 

76.00 

78.00 

79.45 

80.60 

Table  2.--0utput  from  the  sample  problem  (Cont.) 


HEAD 

IN  LAYER 

3  AT  END  OF 

TIME  STEP  1 

IN  STRESS  PERIOD 

1 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

1 

1.795 

24.26 

43.27 

58.61 

71.24 

81.98 

91.41 

99.56 

106.5 

112.2 

117.0 

120.9 

123 . 9 

126.0 

127.1 

2 

1.758 

23.77 

42.37 

57.33 

69.58 

80.00 

89.61 

97 . 92 

104 . 9 

110.5 

115.2 

119.2 

122.3 

124.5 

125.7 

3 

1.685 

22.79 

40.59 

54.79 

66.13 

75.21 

85.92 

94 .71 

101.7 

107.1 

111.5 

115.7 

119.2 

121.7 

122.9 

1.573 

21.29 

37 . 91 

51.09 

60.78 

62.63 

80.35 

90.22 

97.14 

101.8 

104.1 

110.0 

114.5 

117.5 

119.0 

5 

1.412 

19.14 

34.24 

46.68 

57.04 

65.74 

76.48 

85.24 

91.62 

94 . 12 

77.41 

100.6 

108.2 

112.1 

113.9 

6 

1.174 

15.96 

28.86 

40.28 

50.71 

60.62 

70.65 

79.33 

85.96 

90.07 

90.55 

88.50 

101.1 

105.9 

108.0 

7 

0.8258 

11.19 

20.76 

30.85 

41.06 

51.51 

62.64 

72.18 

79.46 

84.42 

87.94 

90.72 

95.89 

99.37 

101.4 

8 

0.4326 

5.126 

10.18 

19.26 

29.18 

39.83 

53.39 

64.03 

72.07 

76.91 

81.53 

84 . 64 

88.83 

91.39 

93.90 

9 

0.7533 

10.21 

18.80 

27.82 

36.64 

45.04 

52.75 

56.99 

64 . 98 

67.60 

73.77 

75.26 

80.67 

81.59 

86.19 

10 

1.037 

14.10 

25.26 

34.81 

42.96 

4  9.62 

54.51 

57.40 

62.57 

65.40 

70.01 

72.29 

76.35 

78.11 

81.39 

11 

1.222 

16.56 

29.33 

39.43 

47.24 

52.75 

55.49 

54 . 97 

60.12 

60.90 

66.29 

67.02 

72.08 

72.56 

77.34 

12 

1.338 

18.11 

31 . 92 

42.49 

50.28 

55.43 

57.90 

58.33 

61.56 

63.04 

66.75 

68.36 

71.93 

73.32 

76.45 

13 

1.412 

19.09 

33.59 

44.56 

52.48 

57.55 

59.59 

58.35 

62.44 

62.4  9 

67.07 

67.30 

71.76 

71.86 

76.20 

14 

1.456 

19.68 

34.62 

45.90 

54.08 

59.58 

62.71 

64.20 

66.82 

68.48 

71.22 

72.87 

75.43 

76.73 

78.67 

15 

1.477 

19.96 

35 . 12 

46.58 

54.95 

60.76 

64.54 

67.06 

69.48 

71.56 

73.83 

75.77 

77.77 

79.22 

80.38 

Table  2.~Output  from  the  sample  problem  (Cont.) 


VOLUMETRIC  BUDGET  FOR  ENTIRE  MODEL  AT  END  OF  TIME  STEP  1  IN  STRESS  PERIOD  1 


CUMULATIVE  VOLUMES  L**3 


RATES  FOR  THIS  TIME  STEP  L**3/T 


K> 

On 


IN: 


IN: 


STORAGE 
CONSTANT  HEAD 
WELLS 
DRAINS 
RECHARGE 


0.00000 
0 .00000 
0.00000 
0.00000 
0 . 13608E+08 


STORAGE  = 
CONSTANT  HEAD  = 
WELLS  = 
DRAINS  = 
RECHARGE  = 


0.00000 
0 . 00000 
0.00000 
0.00000 
157.50 


TOTAL  IN  = 


0 . 13608E+08 


TOTAL  IN  =  157.50 


OUT: 


OUT: 


STORAGE 
CONSTANT  HEAD 
WELLS 
DRAINS 
RECHARGE 


0.00000 
0 . 43291E+07 
0 . 64800E+07 
0.27985E+07 
0.00000 


STORAGE  = 
CONSTANT  HEAD  = 
WELLS  = 
DRAINS  = 
RECHARGE  = 


0 . 00000 
50.105 
75.000 
32.390 
0.00000 


TOTAL  OUT  = 


0 . 13608E+08 


TOTAL  OUT  =  157.50 


IN  -  OUT  = 


406.00 


IN  -  OUT  =  0 . 4 6997E-02 


PERCENT  DISCREPANCY  = 


0.00 


PERCENT  DISCREPANCY  = 


0.00 


TIME  SUMMARY  AT  END  OF  TIME  STEP  1  IN  STRESS  PERIOD  1 

SECONDS  MINUTES  HOURS  DAYS  YEARS 


TIME  STEP  LENGTH 
STRESS  PERIOD  TIME 
TOTAL  SIMULATION  TIME 


86400.0 

1440.00 

24.0000 

1.00000 

0.273785E-02 

86400.0 

1440.00 

24.0000 

1.00000 

0 .273785E-02 

86400.0 

1440.00 

24.0000 

1.00000 

0.273785E-02 

MODULE  DOCUMENTATION 


The  modules  (subroutines)  for  the  GFD  Package  closely  follow  the  design  of  the  BCF  Package 
modules.  GFD  has  one  less  module  becuase  two  modules  were  combined  into  one.  Because  of  the 
similarity  of  design,  including  in  most  cases  identical  comments,  flowcharts  and  detailed  narratives 
are  not  included  here.  Instead,  a  brief  narrative  describes  differences  between  the  corresponding 
GFD  and  BCF  Package  modules.  A  listing  of  each  module  and  a  list  of  variables  are  provided. 

Statements  must  be  added  to  the  MAIN  program  to  call  the  GFD  Package  modules.  There  are  four 
calls  to  GFD  Package  modules  —  GFD1AL,  GFD1RP,  GFD1FM,  and  GFD1BD.  Add  the  GFD 
call  statements  immediately  following  the  corresponding  call  statements  for  the  BCF  Package. 
For  example,  the  GFD1AL  call  should  immediately  follow  the  BCF1 AL  call.  The  calling 
statements  are: 

IF (IUNIT (14) .GT. 0)  CALL  GFDlAL ( ISUM, LENX, LCSC1 , LCCDTR, LCCDTC, 

1  LCBOT , LCTOP , LCSC2 , IUNIT(14) , ISS, NCOL, NROW, NLAY, TOUT, IGFDCB) 

IF (IUNIT (14) .GT. 0)  CALL  GFD1RP (X (LCIBOU) , X (LCHNEW) , X (LCSC1 ) , 

1  X (LCCDTR) , X (LCCDTC) ,X(LCCR) ,X(LCCC) ,X(LCCV) , X (LCDELR) , 

2  X (LCDELC) , X (LCBOT) ,X (LCTOP) , X(LCSC2) , 

3  IUNIT (14) , ISS, NCOL, NROW, NLAY, NODES, IOUT) 

IF (IUNIT (14) .GT.O)  CALL  GFD1FM (X (LCHCOF) , X (LCRHS) , X (LCHOLD) , 

1  X(LCSCl) , X (LCHNEW) , X (LCIBOU) ,X(LCCR) ,X(LCCC) ,X(LCCV) , 

2  X ( LCCDTR) , X ( LCCDTC ) , X { LCBOT ) , X ( LCTOP ) , X ( LCSC2 ) , 

3  DELT, ISS, KKITER,KKSTP,KKPER, NCOL, NROW, NLAY, IOUT) 

IF (IUNIT (14) .GT.O)  CALL  GFD 1BD (VBNM, VBVL, MSUM, X (LCHNEW) , 

1  X (LCIBOU) ,X (LCHOLD) ,X(LCSC1) ,X(LCCR) ,X(LCCC) ,X(LCCV) , 

2  X( LCTOP) , X (LCSC2 ) , DELT, ISS ,  NCOL, NROW, NLAY, KKSTP , KKPER, 

3  IGFDCB , ICBCFL, X ( LCBUFF ) , IOUT ) 

The  above  statements  call  the  GFD  Package  whenever  the  14th  element  of  the  IUNIT  array  is 
greater  than  0.  Another  IUNIT  element  can  be  substituted  if  desired.  To  do  this,  replace 
IUNIT(14)  in  the  “IF”  statements  and  in  the  calling  arguments. 

Modules  GFDlAL  and  GFD1RP  are  very  similar  to  the  corresponding  BCF  Package  modules. 
The  differences  are  in  the  specific  arrays  being  allocated  and  read.  For  example,  an  anisotropy 
array  is  not  needed  in  GFD  because  anisotropy  is  incorporated  in  conductance  terms;  for 
unconfined  and  convertible  layers  conductance  divided  by  thickness  (CDTR  and  CDTC  arrays)  are 
used  by  GFD  instead  of  hydraulic  conductivity  (HY  array). 

Module  GFD1FM  is  essentially  identical  to  the  corrsponding  BCF  Package  module.  The 
changes  all  result  from  a  change  in  a  submodule,  SGFD1H,  which  calculates  conductance  for 
unconfined  and  convertible  layers.  SGFD1H  requires  different  arguments  than  did  the 
analogous  BCF  Package  module  (SBCF1H),  and  these  arguments  must  be  passed  through 
GFD1FM.  The  statements  in  GFD1FM  that  are  not  involved  with  calling  SGFD1H  are  identical 
to  those  in  BCF1FM. 

Module  GFD1BD  is  identical  to  module  BCF1BD  other  than  the  calls  to  submodule.  Submodules 
SGFD1F  and  SGFD1B  are  called  rather  than  SBCF1F  and  SBCF1B. 
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Module  SGFD1N  is  a  simplified  version  of  SBCF1N.  There  is  less  initialization  work  because 
conductance  and  storage  input  parameters  more  closely  match  the  fundamental  data  used  in  the 
flow  equation. 

Module  SGFD1H  calculates  horizontal  conductance  for  unconfined  and  convertible  layers.  First 
it  calculates  saturated  thickness  using  aquifer  bottom  elevation,  aquifer  top  elevation  if  appropriate, 
and  the  current  head.  Then  it  calculates  conductance  from  saturated  thickness  and  conductance 
divided  by  thickness.  This  module  provides  the  analogous  functions  of  modules  SBCF1H  and 
SBCF1C  in  the  BCF  Package.  SBCF1H  calculates  saturated  thickness  and  transmissivity; 
SBCF1C  calculates  conductance  from  transmissivity  and  cell  dimensions.  SBCF1C  was  separated 
from  SBCF1H  because  the  calculation  of  conductance  from  transmissivity  was  also  required  by 
module  SBCF1N  for  confined  layers.  By  putting  the  transmissivity  calculation  in  a  separate 
module,  duplication  of  this  code  was  avoided.  The  GFD  Package  does  not  require  the  calculation 
of  conductance  by  module  SGFD1N  for  confined  layers.  SGFD1H  is  the  only  module  that  requires 
conductance  calculations.  Thus,  there  is  no  reason  to  have  separate  modules  for  calculating 
saturated  thickness  and  conductance.  The  steps  and  structure  in  SGFD1H  closely  match  those  of 
routines  SBCF1H  and  S BCF  1C  combined. 

Modules  SGFD1B  and  SGFD1F  are  identical  to  SBCF1B  and  SBCF1F,  except  that  variable 
1GFDCB  replaces  variable  EBCFCB. 

Memoiy  requirements  for  GFD  data  are  similar  to  BCF  memory  requirements.  Most  data  in  the 
model  are  stored  in  the  X  array  (McDonald  and  Harbaugh,  1988,  p.  3-22  and  Appendix  B-l),  and 
the  dimensioned  size  of  the  X  array  must  be  large  enough  to  hold  the  data  used  by  all  packages. 
The  X-array  requirements  for  the  GFD  Package  can  be  calculated  by  adding  the  space  required  for 
each  GFD  option. 

A.  If  a  transient  simulation  (ISS  is  0),  add  NCOL*NROW*NLAY. 

B.  For  each  layer  where  LAYCON  is  one  or  three,  add  3*NCOL*NROW. 

C.  For  each  layer  where  LAYCON  is  two  or  three,  add  NCOL*NROW. 

D.  If  a  transient  simulation  (ISS  is  0),  for  each  layer  where  LAYCON  is  two  or  three,  add 
NCOL*NROW. 
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GFD1AL  Module 


SUBROUTINE  GFD1AL (ISUM, LENX, LCSC1 , LCCDTR, LCCDTC, LCBOT, 

1  LCTOP , LCSC2 , IN, ISS,NCOL,NROW,NLAY, IOUT, IGFDCB) 

C 

C - VERSION  1304  19SEP1989  GFD1AL 

C 

Q  ****************************************************************** 

C  ALLOCATE  ARRAY  STORAGE  FOR  GENERAL  FINITE-DIFFERENCE  FLOW  PACKAGE 

Q  ****************************************************************** 

c 

C  SPECIFICATIONS: 

C  - 

COMMON  /FLWCOM/ LAYCON (80) 

c  - 

C 

Cl - IDENTIFY  PACKAGE 

WRITE (IOUT, 1) IN 

1  FORMAT ( 1H0 , ' GFD1  —  GENERAL  FINITE-DIFFERENCE  FLOW  PACKAGE,  ' , 

1  ' VERSION  1,  9/19/89  INPUT  READ  FROM  UNIT',13) 

C 

C2 - READ  AND  PRINT  ISS  (STEADY-STATE  FLAG)  AND  IGFDCB  (FLAG  FOR 

C2 - PRINTING  OR  UNIT#  FOR  RECORDING  CELL-BY-CELL  FLOW  TERMS) 

READ (IN, 2)  ISS, IGFDCB 

2  FORMAT (2110) 

IF(ISS.EQ.O)  WRITE (IOUT, 3) 

3  FORMAT (IX, ' TRANSIENT  SIMULATION') 

IF(ISS.NE.O)  WRITE (IOUT, 4) 

4  FORMAT (IX, 'STEADY-STATE  SIMULATION' ) 

IF (IGFDCB. GT. 0)  WRITE ( IOUT, 9 )  IGFDCB 

9  FORMAT (IX, 'CELL  BUDGET  WILL  BE  SAVED  ON  UNIT',13) 

IF (IGFDCB.LT. 0)  WRITE (IOUT, 88) 

88  FORMAT (IX, 'CONSTANT  HEAD  CELL-BY-CELL  FLOWS  WILL  BE  PRINTED') 

C 

C3 - READ  TYPE  CODE  FOR  EACH  LAYER  AND  COUNT  TOPS  AND  BOTTOMS 

IF (NLAY. LE . 80)  GO  TO  50 
WRITE (IOUT, 11) 

11  FORMAT (1 HO, ' YOU  HAVE  SPECIFIED  MORE  THAN  80  MODEL  LAYERS' /IX, 

1  'SPACE  IS  RESERVED  FOR  A  MAXIMUM  OF  80  LAYERS  IN  ARRAY  LAYCON') 
STOP 
C 

C3A — - — READ  LAYER  TYPE  CODES. 

50  READ (IN, 51)  (LAYCON (I) , 1=1, NLAY) 

51  FORMAT ( 4012 ) 

C  BOTTOM  IS  READ  FOR  TYPES  1,3  TOP  IS  READ  FOR  TYPES  2,3 

WRITE (IOUT, 52) 

52  FORMAT (IX, 5X, ' LAYER  AQUIFER  TYPE' ,/ IX, 5X, 1 9 ('-') ) 

C 

C3B - INITIALIZE  TOP  AND  BOTTOM  COUNTERS. 

NBOT=0 

NTOP=0 

C 

C3C - PRINT  LAYER  TYPE  AND  COUNT  TOPS  AND  BOTTOMS  NEEDED. 

DO  100  1=1, NLAY 
C 

C3C1 - PRINT  LAYER  NUMBER  AND  LAYER  TYPE  CODE. 

L=LAYCON ( I ) 

WRITE (IOUT, 7)  I , L 

7  FORMAT (IX, 19, 110) 


29 


c 

C3C2 - ONLY  THE  TOP  LAYER  CAN  BE  UNCONFINED (LAYCON=l) . 

IF ( L . NE . 1  .OR.  I.EQ.l)  GO  TO  70 
WRITE ( IOUT ,  8) 

8  FORMAT (1H0, 'AQUIFER  TYPE  1  IS  ONLY  ALLOWED  IN  TOP  LAYER') 
STOP 
C 

C3C3 - LAYER  TYPES  1  AND  3  NEED  A  BOTTOM.  ADD  1  TO  KB. 

70  IF ( L . EQ . 1  .OR.  L.EQ.3)  NBOT=NBOT+l 
C 

C3C4 - LAYER  TYPES  2  AND  3  NEED  A  TOP.  ADD  1  TO  KT . 

IF (L . EQ . 2  .OR.  L.EQ.3)  NTOP=NTOP+l 
100  CONTINUE 


C 

C 

C 

C4 - COMPUTE  DIMENSIONS  FOR  ARRAYS. 

NRC=NROW*NCOL 
IS IZ=NRC*NLAY 
C 

C5 - ALLOCATE  SPACE  FOR  ARRAYS.  IF  RUN  IS  TRANSIENT ( ISS=0 ) 

C5 - THEN  SPACE  MUST  BE  ALLOCATED  FOR  STORAGE. 

ISOLD=ISUM 

LCSC1=ISUM 

IF(ISS.EQ.O)  ISUM-ISUM+ISIZ 
*  .  LCSC2  =  ISUM 


' IF(ISS.EQ.O)  I SUM- 1 SUM+NRC  *NTOP 
LCBOT— ISUM 
ISUM= I SUM+NRC  *NBOT 
LCCDTR=ISUM 
ISUM= I SUM+NRC  *NBOT 
LCCDTC=ISUM 
ISUM= I SUM+NRC  *NBOT 
LCTOP=ISUM 
I  S'UM=  I  SUM+NRC  *NTOP 
C 

C6 - PRINT  THE  AMOUNT  OF  SPACE  USED  BY  THE  GFD  PACKAGE. 

I SP= I SUM- I SOLD 
WRITE (IOUT, 101)  ISP 

101  FORMAT (IX, 18 , '  ELEMENTS  IN  X  ARRAY  ARE  USED  BY  GFD') 
ISUM1=ISUM-1 

WRITE (IOUT, 102)  ISUMl , LENX 

102  FORMAT (IX, 18, '  ELEMENTS  OF  X  ARRAY  USED  OUT  OF',18) 

IF (ISUMl .GT. LENX)  WRITE (IOUT, 103) 

103  FORMAT (IX, '  ***X  ARRAY  MUST  BE  DIMENSIONED  LARGER***') 

C 

C7 - RETURN 

RETURN 

END 
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List  of  Variables  for  Module  GFD1AL 


Variable 


I 

IGFDCB 


IN 

IOUT 

ISIZ 

ISOLD 

ISP 

ISS 


I  SUM 


ISUM1 

L 

LAYCON 


LCBOT 

LCCDTC 

LCCDTR 

LCSC1 

LCSC2 

LCTOP 

LENX 


NBOT 

NCOL 

NLAY 

NRC 

NROW 

NTOP 


Range 


Definition 


Module 

Package 


Package 

Global 

Module 

Module 

Module 

Package 


Global 


Module 

Module 

Package 


Package 

Package 

Package 

Package 

Package 

Package 

Global 


Module 

Global 

Global 

Module 

Global 

Module 


Index 

Flag  and  a  unit  number. 

>  0,  unit  number  on  which  cell-by-cell  flow  terms 
will  be  recorded  whenever  ICBCFL  is  set. 

=  0f  cell-by-cell  flow  terms  will  not  be  printed  or 
recorded. 

<  0,  flow  from  each  constant-head  cell  will  be  printed 
whenever  ICBCFL  is  set. 

Primary  unit  number  from  which  input  for  this  package  will 
be  read. 

Primary  unit  number  for  all  printed  output. 

Number  of  cells  in  the  grid. 

Value  of  ISUM  upon  entry  to  this  module. 

Number  of  elements  in  the  X  array  allocated  by  this  module. 
Steady-state  flag. 

=  0,  simulation  is  transient 
not  0,  simulation  is  steady  state. 

Element  number  of  the  lowest  element  in  the  X  array  that  has 
not  yet  been  allocated.  When  space  is  allocated  in  the 
X  array,  the  size  of  the  allocation  is  added  to  ISUM. 
ISUM-1 

Temporary  equivalent  of  LAYCON (I). 

DIMENSION  (80)  Layer-type  code: 

0  —  Layer  is  strictly  confined. 

1  —  Layer  is  strictly  unconfined. 

2  —  Layer  is  convertible  between  confined  and  unconfined 

(saturated  thickness  is  constant) . 

3  --  Layer  is  convertible  between  confined  and  unconfined 

(saturated  thickness  varies) . 

Location  in  the  X  array  of  the  first  element  of  array  BOT . 

Location  in  the  X  array  of  the  first  element  of  array  CDTC . 

Location  in  the  X  array  of  the  first  element  of  array  CDTR. 

Location  in  the  X  array  of  the  first  element  of  array  SCI. 

Location  in  the  X  array  of  the  first  element  of  array  SC2 . 

Location  in  the  X  array  of  the  first  element  of  array  TOP. 

Number  of  elements  in  the  X  array.  This  should  always  be 
equal  to  the  dimension  of  X  specified  in  the  MAIN 
program. 

Counter  for  the  number  of  layers  that  need  elevation  of  the 
bottom  and  conductance  divided  by  thickness  arrays . 
Number  of  columns  in  the  grid. 

Number  of  layers  in  the  grid. 

Number  of  cells  in  a  layer. 

Number  of  rows  in  the  grid. 

Counter  for  the  number  of  layers  that  need  elevation  of  the 
top  and  (if  transient)  secondary  storage  capacity 
arrays . 
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GFD1RP  Module 


SUBROUTINE  GFD1RP ( IBOUND , HNEW, SCI , CDTR, CDTC, CR, CC, CV, DELR, DELC, 

1  BOT, TOP , SC2 , IN, ISS , NCOL, NROW, NLAY, NODES , IOUT) 

C 

C - VERSION  1406  19SEP198  9  GFD1RP 

C 

0  ****************************************************************** 

C  READ  AND  INITIALIZE  DATA  FOR  GENERAL  FLOW  PACKAGE 

0  ****************************************************************** 

C 

C  SPECIFICATIONS: 

c  - 

CHARACTER* 4  ANAME 
DOUBLE  PRECISION  HNEW 
C 

DIMENSION  HNEW (NODES) , SCI (NODES) , CR (NODES) , CC (NODES) , CV (NODES) , 

1  CDTR (NODES) , CDTC (NODES) , ANAME (6, 11) , DELR (NCOL) , 

2  DELC ( NROW ) , BOT (NODES) , TOP (NODES) , SC2 (NODES) , 

3  IBOUND (NODES) 

C 

COMMON  /FLWCOM/ LAYCON (80) 

C 

DATA  ANAME (1,1), ANAME (2,1), ANAME (3,1), ANAME (4,1), ANAME (5,1), 

1  ANAME (6,1)  / ' PRIM' , ' ARY  ' , ' STOR' , ' AGE  CAPA' ,' CITY' / 

DATA  ANAME (1,2) , ANAME (2, 2) , ANAME (3, 2) ,ANAME(4,2) ,ANAME(5,2) , 

1  ANAME (6,2)  /'  CO ' , ' NDUC '  ,  ' T ANC ' , ' E  AL ' , ' ONG  ', 'ROWS'/ 

DATA  ANAME (1,3) ,ANAME(2,3) , ANAME (3, 3) ,ANAME(4,3) ,ANAME(5,3) , 

1  ANAME (6,3)  /'  C' , ' OND/' , ' THIC' , ' K  AL' , ' ONG  ', 'ROWS'/ 

DATA  ANAME (1, 4) , ANAME (2, 4) , ANAME (3,4) , ANAME (4, 4) , ANAME (5, 4) , 

1  ANAME (6, 4)  /'  ' , ' VERT' , ' I CAL' , '  CON' , ' DUCT' , ' ANCE' / 

DATA  ANAME (1, 5) , ANAME (2,5) , ANAME (3, 5) , ANAME (4, 5) , ANAME (5, 5) , 

1  ANAME (6, 5)  /'  ','  ','  ','  ','  BO' , ' TTOM' / 

DATA  ANAME (1,6) , ANAME (2,6) , ANAME (3,6) , ANAME (4,6) , ANAME (5,6), 

1  ANAME (6, 6)  /'  ','  ','  ','  ','  ','  TOP'/ 

DATA  ANAME (1, 7) , ANAME (2,7) , ANAME (3,7) , ANAME (4, 7) , ANAME (5,7), 

1  ANAME (6, 7)  /'  S' , ' EC .  ' , ' STOR' , ' AGE  ',' CAPA' ,' CITY' / 

DATA  ANAME (1, 8) , ANAME (2,8) , ANAME (3, 8) , ANAME (4, 8) , ANAME (5, 8) , 

1  ANAME (6,8)  /'  CO' , ' NDUC' , ' TANC' , ' E  AL' , ' ONG  ', 'COLS'/ 

DATA  ANAME (1, 9) , ANAME (2, 9) , ANAME (3, 9) , ANAME (4, 9) , ANAME (5, 9) , 

1  ANAME (6, 9)  /'  ','  ','  ','  ' , ' DELR' / 

DATA  ANAME (1, 10) , ANAME (2,10) , ANAME (3,10) , ANAME (4, 10) , ANAME (5, 10) , 

1  ANAME (6, 10)  /'  ','  ','  ','  ','DELC'/ 

DATA  ANAME (1,11), ANAME (2,11), ANAME (3,11), ANAME (4,11), ANAME (5,11), 

1  ANAME (6, 11)  /'  C' , ' OND/ ' , ' THIC' , ' K  AL','ONG  ', 'COLS'/ 


Cl - CALCULATE  NUMBER  OF  NODES  IN  A  LAYER  AND  READ  DELR, DELC 

NI J=NCOL*NROW 
C 

CALL  U1DREL (DELR, ANAME (1, 9) ,NCOL, IN, IOUT) 

CALL  U1DREL (DELC, ANAME ( 1, 10) ,NROW, IN, IOUT) 

C 

02-- - READ  ALL  PARAMETERS  FOR  EACH  LAYER 

KT=0 

KB=0 

DO  200  K=1 , NLAY 
KK=K 
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c 

C2A - FIND  ADDRESS  OF  EACH  LAYER  IN  THREE  DIMENSION  ARRAYS. 

IF (LAYCON(K) .EQ.l  .OR.  LAYCON (K) . EQ . 3)  KB=KB+1 
IF (LAYCON(K) .EQ.2  .OR.  LAYCON (K) . EQ . 3 )  KT=KT+1 
LOC=l+ (K-l ) *NIJ 
LOCB=l+ (KB-1) *NI  J 
LOCT=l+ (KT-1 ) *NI  J 


C 

C2B - READ  PRIMARY  STORAGE  CAPACITY  IF  TRANSIENT  SIMULATION 

IF ( ISS . EQ . 0 ) CALL  U2DREL ( SCI (LOC) , ANAME (1,1) , NROW, NCOL, KK, IN, IOUT) 
C 

C2C - READ  CONDUCTANCE  IF  LAYER  TYPE  IS  0  OR  2 

IF (LAYCON (K) .EQ. 3  .OR.  LAYCON (K) . EQ . 1)  GO  TO  100 


CALL  U2DREL(CR(LOC) , ANAME < 1,2)  , NROW, NCOL, KK, IN,  IOUT) 
CALL  U2DREL (CC (LOC) , ANAME (1, 8) , NROW, NCOL, KK, IN, IOUT) 
GO  TO  110 


C 

C2D - READ  SPECIFIC  CONDUCTANCE  AND  BOTTOM  ELEVATION (BOT) 

C2D - IF  LAYER  TYPE  IS  1  OR  3 


100  CALL  U2DREL (CDTR(LOCB) , ANAME (1, 3) , NROW, NCOL, KK, IN, IOUT) 
CALL  U2DREL (CDTC (LOCB) , ANAME (1, 11) , NROW, NCOL, KK, IN, IOUT) 
CALL  U2DREL(BOT (LOCB) , ANAME (1,5) , NROW, NCOL, KK, IN, IOUT) 


C 

C2E - READ  VERTICAL  CONDUCTANCE  IF  NOT  BOTTOM  LAYER 

110  IF (K.EQ.NLAY)  GO  TO  120 

CALL  U2DREL (CV (LOC) , ANAME ( 1 , 4) , NROW, NCOL, KK, IN, IOUT) 

C 

C2F - READ  SECONDARY  STORAGE  CAPACITY  IF  TRANSIENT  SIMULATION 

C2F - AND  LAYER  TYPE  IS  2  OR  3 

120  IF (LAYCON (K) .NE. 3  .AND.  LAYCON (K) .NE . 2 )  GO  TO  200 

IF ( ISS . EQ . 0 ) CALL  U2DREL(SC2 (LOCT) ,ANAME(1,7) , NROW, NCOL, KK, IN, IOUT) 
C 

C2G - READ  TOP  ELEVATION (TOP )  IF  LAYER  TYPE  IS  2  OR  3 

CALL  U2DREL(TOP (LOCT) , ANAME (1, 6) , NROW, NCOL, KK, IN, IOUT) 

200  CONTINUE 
C 

C3 - INITIALIZE  AND  CHECK  GFD  DATA 

CALL  SGFDlN (HNEW, IBOUND, CR, CC, CV, CDTR, CDTC, 

1  NCOL, NROW, NLAY, IOUT) 

C 

C4 - RETURN 

RETURN 

END 
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List  of  Variables  for  Module  GFD1RP 
Variable  Range  Definition 


ANAME  Module  DIMENSION  (6,11),  Labels  for  printout  of  input  arrays. 

BOT  Package  DIMENSION  (NODES) ,  Elevation  of  the  aquifer  bottom. 

Although  BOT  is  dimensioned  to  the  size  of  the  grid, 
space  exists  only  for  cells  for  which  saturated 
thickness  is  calculated. 

CC  Global  DIMENSION  (NODES),  Conductance  in  the  column  direction. 

CDTC  Package  DIMENSION  (NODES) ,  Conductance  divided  by  thickness  in  the 

column  direction.  Although  CDTC  is  dimensioned  to  the 
size  of  the  grid,  space  exists  only  for  cells  for  which 
saturated  thickness  is  calculated. 

CDTR  Package  DIMENSION  (NODES) ,  Conductance  divided  by  thickness  in  the 

row  direction.  Although  CDTR  is  dimensioned  to  the 
size  of  the  grid,  space  exists  only  for  cells  for  which 
saturated  thickness  is  calculated. 

CR  Global  DIMENSION  (NODES),  Conductance  in  the  row  direction. 

CV  Global  DIMENSION  (NODES),  Conductance  in  the  vertical  direction. 

DELC  Global  DIMENSION  (NROW) ,  Cell  dimension  in  the  column  direction. 

DELC(I)  contains  the  width  of  row  I. 

DELR  Global  DIMENSION  (NCOL) ,  Cell  dimension  in  the  row  direction. 

DELR(J)  contains  the  width  of  column  J. 

HNEW  Global  DIMENSION  (NODES) ,  Most  recent  estimate  of  head  in  each 

cell . 

I BOUND  Global  DIMENSION  (NODES),  Status  of  each  cell: 

<  0,  constant-head  cell 
=  0,  no-flow  cell 
>  0,  variable-head  cell 

IN  Package  Primary  unit  number  from  which  input  for  this  package  will 

be  read. 

IOUT  Global  Primary  unit  number  for  all  printed  output. 

ISS  Package  Steady-state  flag. 

=  0,  simulation  is  transient 

not  0,  simulation  is  steady  state. 

K  Module  Index  for  layers . 

KB  Module  Index  for  layer  within  BOT,  CDTR,  and  CDTC  arrays. 

KK  Module  Temporary  variable  set  equal  to  K.  KK  is  used  as  an  actual 

argument  in  subroutine  calls  to  avoid  using  the  DO  loop 
variable  K  as  an  argument,  which  causes  problems  with 
some  compilers. 

KT  Module  Index  for  layer  within  TOP  and  (if  transient)  SC2  arrays. 

LAYCON  Package  DIMENSION  (80)  Layer-type  code: 

0  —  Layer  is  strictly  confined. 

1  --  Layer  is  strictly  unconfined. 

2  —  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  is  constant) . 

3  --  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  varies) . 


LOC 

Module 

Pointer  used  to  point  to  different 

layers 

in 

conductance 

arrays . 

LOCB 

Module 

Pointer  used  to  point  to  different 

layers 

in 

BOT,  CDTR,  and 

CDTC  arrays . 

LOCT 

Module 

Pointer  used  to  point  to  different 

layers 

in 

TOP  and  SC2 

arrays . 

NCOL 

Global 

Number  of  columns  in  the  grid. 

NIJ 

Module 

Number  of  cells  in  a  layer. 

NLAY 

Global 

Number  of  layers  in  the  grid. 
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NODES 

NROW 

SCI 

SC2 


TOP 


Global 

Global 

Package 

Package 


Package 


Number  of  cells  (nodes)  in  the  finite-difference  grid. 
Number  of  rows  in  the  grid. 

DIMENSION  (NODES) ,  Primary  storage  capacity. 

DIMENSION  (NODES) ,  Secondary  storage  capacity. 

Although  SC2  is  dimensioned  to  the  size  of  the  grid, 
space  exists  only  for  cells  that  can  convert  between 
confined  and  unconfined. 

DIMENSION  (NODES),  Elevation  of  the  aquifer  top. 

Although  TOP  is  dimensioned  to  the  size  of  the  grid, 
space  exists  only  for  cells  that  can  convert  between 
confined  and  unconfined. 
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GFD1FM  Module 


SUBROUTINE  GFD1FM (HCOF, RHS, HOLD, SC1,HNEW, IBOUND, CR, CC,  CV, 

1  CDTR, CDTC, BOT, TOP , SC2 , DELT, ISS, KITER,  KSTP,  KPER, 

2  NCOL, NROW, NLAY, IOUT) 

C - VERSION  0912  1 9SEP198  9  GFD1FM 

C 

Q  ********************************************************************** 

C  ADD  LEAKAGE  CORRECTION  AND  STORAGE  TO  HCOF  AND  RHS,  AND  CALCULATE 

C  CONDUCTANCE  AS  REQUIRED 

0  ******************************************************************** 

C 

C  SPECIFICATIONS: 

c  - 

DOUBLE  PRECISION  KNEW 
C 

DIMENSION  HCOF (NCOL, NROW, NLAY) , RHS (NCOL, NROW, NLAY) , 

1  HOLD (NCOL, NROW, NLAY) , SCI (NCOL, NROW, NLAY) , HNEW (NCOL, NROW, NLAY) , 

2  IBOUND (NCOL, NROW, NLAY) , CR (NCOL, NROW, NLAY) , 

3  CC (NCOL, NROW, NLAY) , CV (NCOL, NROW, NLAY) , CDTR (NCOL, NROW, NLAY) , 

4  CDTC (NCOL, NROW, NLAY) , BOT (NCOL, NROW, NLAY) , TOP (NCOL, NROW, NLAY) , 

5  SC2 (NCOL, NROW, NLAY) 

C 

COMMON  / FLWCOM/ LAYCON (80) 

c  - 

KB=0 

KT=0 

C 

Cl - FOR  EACH  LAYER:  IF  T  VARIES  CALCULATE  HORIZONTAL  CONDUCTANCES 

DO  100  K=1 , NLAY 
KK=K 

IF (LAYCON (K) . EQ . 3  .OR.  LAYCON (K) . EQ . 2 )  KT=KT+1 
C 

CIA - IF  LAYER  TYPE  IS  NOT  1  OR  3  THEN  SKIP  THIS  LAYER. 

IF (LAYCON (K) . NE . 3  .AND.  LAYCON (K) . NE . 1 )  GO  TO  100 
KB=KB+1 
C 

C1B - FOR  LAYER  TYPES  1  &  3  CALL  SGFDC1  TO  CALCULATE 

C1B - HORIZONTAL  CONDUCTANCES. 

CALL  SGFD1H (HNEW, IBOUND, CR, CC, CV, CDTR, CDTC, BOT, TOP, 

1  KK, KB, KT,KITER, KSTP, KPER, NCOL, NROW, NLAY, IOUT) 

100  CONTINUE 
C 

C2 - IF  THE  SIMULATION  IS  TRANSIENT  ADD  STORAGE  TO  HCOF  AND  RHS 

IF (ISS .NE. 0)  GO  TO  201 

TLED=1./DELT 

KT=0 

DO  200  K=1 , NLAY 
C 

C3 - SEE  IF  THIS  LAYER  IS  CONVERTIBLE  OR  NON-CONVERTIBLE. 

IF (LAYCON (K) . EQ . 3  .OR.  LAYCON (K) . EQ . 2 )  GO  TO  150 

C4 - NON-CONVERTIBLE  LAYER,  SO  USE  PRIMARY  STORAGE 

DO  140  1=1, NROW 
DO  140  J=1 , NCOL 

IF (IBOUND (J, I,K) .LE. 0)  GO  TO  140 

RHO=SCl (J, I , K) *TLED 

HCOF ( J, I , K) =HCOF ( J, I, K) -RHO 

RHS ( J, I, K) =RHS (J, I , K) -RHO*HOLD ( J, I, K) 

140  CONTINUE 
GO  TO  200 
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c 

C5 - A  CONVERTIBLE  LAYER,  SO  CHECK  OLD  AND  NEW  HEADS  TO  DETERMINE 

C5 - WHEN  TO  USE  PRIMARY  AND  SECONDARY  STORAGE 

150  KT=KT+1 

DO  180  1=1 , NROW 
DO  180  J=1 , NCOL 
C 

C5A - IF  THE  CELL  IS  EXTERNAL  THEN  SKIP  IT. 

IF (IBOUND (J, I,K) .LE. 0)  GO  TO  180 
TP=TOP ( J, I, KT) 

RH02=SC2 (J, I,KT) *TLED 
RH01=SC1 (J, I,K) *TLED 
C 

C5B - FIND  STORAGE  FACTOR  AT  START  OF  TIME  STEP. 

SOLD=RH02 

IF (HOLD (J, I,K) .GT.TP)  SOLD=RH01 
C 

C5C - FIND  STORAGE  FACTOR  AT  END  OF  TIME  STEP. 

HTMP=HNEW ( J, I , K) 

SNEW=RH02 

IF (HTMP .GT .TP)  SNEW=RH01 
C 

C5D - ADD  STORAGE  TERMS  TO  RHS  AND  HCOF . 

HCOF ( J, I, K) =HCOF ( J, I, K) -SNEW 

RHS (J, I,K)=RHS (J, I, K)  -  SOLD* (HOLD (J, I, K) -TP)  -  SNEW*TP 
C 

180  CONTINUE 
C 

200  CONTINUE 
C 

C6 - FOR  EACH  LAYER  DETERMINE  IF  CORRECTION  TERMS  ARE  NEEDED  FOR 

C6 - FLOW  DOWN  INTO  PARTIALLY  SATURATED  LAYERS. 

201  KT=0 

DO  300  K=1 , NLAY 
C 

C7 - SEE  IF  CORRECTION  IS  NEEDED  FOR  LEAKAGE  FROM  ABOVE. 

IF (LAYCON (K) . NE.3  .AND.  LAYCON (K) . NE . 2 )  GO  TO  250 
KT=KT+1 

IF(K.EQ.l)  GO  TO  250 
C 

C7A - FOR  EACH  CELL  MAKE  THE  CORRECTION  IF  NEEDED. 

DO  220  1=1 f NROW 
DO  220  J=1 , NCOL 
C 

C7B - IF  THE  CELL  IS  EXTERNAL ( IBOUND<=0 )  THEN  SKIP  IT. 

IF (IBOUND (J, I,K) .LE.0)  GO  TO  220 
HTMP=HNEW ( J, I, K) 

C 

C7C - IF  HEAD  IS  ABOVE  TOP  THEN  CORRECTION  NOT  NEEDED 

IF (HTMP . GE . TOP ( J, I, KT) )  GO  TO  220 
C 

C7D - WITH  HEAD  BELOW  TOP  ADD  CORRECTION  TERMS  TO  RHS  AND  HCOF. 

RHS (Jf I,K)=RHS (J, I,K)  +  CV ( J, I, K-l ) *TOP ( J, I, KT) 

HCOF (J, I,K)=HCOF(J, I,K)  +  CV(J,I,K-1) 

220  CONTINUE 
C 


37 


C8 - SEE  IF  THIS  LAYER  MAY  NEED  CORRECTION  FOR  LEAKAGE  TO  BELOW. 

250  IF (K.EQ.NLAY)  GO  TO  300 

IF (LAYCON (K+l ) . NE . 3  .AND.  LAYCON (K+l ) . NE . 2 )  GO  TO  300 
KTT=KT+1 
C 

C8A - FOR  EACH  CELL  MAKE  THE  CORRECTION  IF  NEEDED. 

DO  280  1=1 , NROW 
DO  280  J=l,NCOL 
C 

C8B - IF  CELL  IS  EXTERNAL  (IBOUND<=0)  THEN  SKIP  IT. 

IF ( IBOUND ( J, I, K) . LE . 0 )  GO  TO  280 
C 

C8C - -IF  HEAD  IN  THE  LOWER  CELL  IS  LESS  THAN  TOP  ADD  CORRECTION 

C8C - TERM  TO  RHS  . 

HTMP=HNEW ( J, I, K+l) 

IF (HTMP . LT . TOP (J, 1, KTT) )  RHS (J, I,K)=RHS (J, I,K) 

1  -  CV(J, I,K) * (TOP (J, I, KTT) -HTMP) 

280  CONTINUE 
300  CONTINUE 
C 

C9 - RETURN 

RETURN 

END 
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List  of  Variables  for  Module  GFD1FM 


Variable  Range 


Definition 


BOT 


CC 

CDTC 


CDTR 


CR 

CV 

DELT 

HCOF 

HNEW 

HOLD 

HTMP 

I 

I BOUND 


IOUT 

ISS 


J 

K 

KB 

KITER 

KK 


KPER 

KSTP 

KT 

KTT 

LAYCON 


Package 


Global 

Package 


Package 


Global 

Global 

GLOBAL 

Global 

Global 

Global 

Module 

Module 

Global 


Global 

Package 


Module 

Module 

Module 

Global 

Module 


Global 

Global 

Module 

Module 

Package 


DIMENSION  (NCOL, NROW, NLAY) ,  Elevation  of  the  aquifer 

bottom.  Although  BOT  is  dimensioned  to  the  size  of  the 
grid,  space  exists  only  for  cells  for  which  saturated 
thickness  is  calculated. 

DIMENSION  (NCOL, NROW, NLAY)  ,  Conductance  in  the  column 
direction . 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  divided  by 

thickness  in  the  column  direction.  Although  CDTC  is 
dimensioned  to  the  size  of  the  grid,  space  exists  only 
for  cells  for  which  saturated  thickness  is  calculated. 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  divided  by 
thickness  in  the  row  direction.  Although  CDTR  is 
dimensioned  to  the  size  of  the  grid,  space  exists  only 
for  cells  for  which  saturated  thickness  is  calculated. 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  row 
direction . 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  vertical 
direction . 

Length  of  the  current  time  step. 

DIMENSION  (NCOL, NROW, NLAY) ,  Coefficient  of  head  in  the 
finite-difference  equation. 

DIMENSION  (NCOL, NROW, NLAY) ,  Most  recent  estimate  of  head  in 
each  cell. 

DIMENSION  (NCOL, NROW, NLAY)  ,  Head  at  the  start  of  the  current 
time  step. 

Temporary  single  precision  equivalent  of  HNEW(J,I,K). 

Index  for  rows . 

DIMENSION  (NCOL, NROW, NLAY) ,  Status  of  each  cell: 

<  0,  constant-head  cell 
=  0,  no-flow  cell 
>  0,  variable-head  cell 

Primary  unit  number  for  all  printed  output. 

Steady-state  flag. 

=  0,  simulation  is  transient 

not  0,  simulation  is  steady  state. 

Index  for  columns . 

Index  for  layers . 

Index  for  layer  within  BOT,  CDTR,  and  CDTC  arrays. 

Iteration  counter.  KITER=1  at  the  start  of  each  time  step. 

Temporary  variable  set  equal  to  K.  KK  is  used  as  an  actual 
argument  in  subroutine  calls  to  avoid  using  the  DO  loop 
variable  K  as  an  argument,  which  causes  problems  with 
some  compilers . 

Stress  period  counter. 

Time  step  counter.  KSTP=1  at  the  start  of  each  stress 
period . 

Index  for  layer  within  TOP  and  (if  transient)  SC2  arrays. 

Index  to  TOP  array  of  layer  immediately  below  layer  K. 

DIMENSION  (80)  Layer-type  code: 

0  --  Layer  is  strictly  confined. 

1  --  Layer  is  strictly  unconfined. 

2  --  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  is  constant) . 

3  —  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  varies) . 
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NCOL  Global  Number  of  columns  in  the  grid. 

NLAY  Global  Number  of  layers  in  the  grid. 

NROW  Global  Number  of  rows  in  the  grid. 

RHO  Module  Storage  capacity  divided  by  time  step  length  for  strictly 

confined  or  strictly  unconfined  layers. 

RHOl  Module  Confined  storage  capacity  divided  by  time  step  length  for 

convertible  layers. 

RH02  Module  Unconfined  storage  capacity  divided  by  time  step  length  for 

convertible  layers. 

RHS  Global  DIMENSION  (NCOL, NROW, NLAY) ,  Right-hand  side  of  finite- 

difference  equation. 

SCI  Package  DIMENSION  (NCOL, NROW, NLAY) ,  Primary  storage  capacity. 

SC2  Package  DIMENSION  (NCOL, NROW, NLAY) ,  Secondary  storage  capacity. 

Although  SC2  is  dimensioned  to  the  size  of  the  grid, 
space  exists  only  for  cells  that  can  convert  between 
confined  and  unconfined. 

SNEW  Module  Storage  capacity  divided  by  time  step  length  at  the  end  of 

the  time  step  for  convertible  layers. 

SOLD  Module  Storage  capacity  divided  by  time  step  length  at  the  start 

of  the  time  step  for  convertible  layers . 

TLED  Module  1/DELT . 

TOP  Package  DIMENSION  (NCOL, NROW, NLAY) ,  Elevation  of  the  aquifer  top. 

Although  TOP  is  dimensioned  to  the  size  of  the  grid, 
space  exists  only  for  cells  that  can  convert  between 
confined  and  unconfined. 

TP  Module  Temporary  equivalent  of  TOP(J,I,KT). 
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GFD1BD  Module 


SUBROUTINE  GFD1BD (VBNM, VBVL, MSUM, HNEW, IBOUND, HOLD, SCI, CR, CC, CV, 

1  TOP, SC2,DELT, ISS , NCOL, NROW, NLAY, KSTP , KPER, IGFDCB, 

2  ICBCFL,BUFF, IOUT) 

C - VERSION  1346  19SEP1989  GFD1BD 

C 

Q  ****************************************************************** 

C  COMPUTE  BUDGET  FLOW  TERMS  FOR  GFD  —  STORAGE,  CONSTANT  HEAD,  AND 

C  FLOW  ACROSS  CELL  WALLS 

0  ****************************************************************** 

C 

C  SPECIFICATIONS: 

C  - 

CHARACTER* 4  VBNM, TEXT 
DOUBLE  PRECISION  HNEW 
C 

DIMENSION  HNEW (NCOL,NROW, NLAY) ,  IBOUND (NCOL, NROW, NLAY) , 

1  HOLD (NCOL, NROW, NLAY) ,  SCI (NCOL, NROW, NLAY) , 

2  CR (NCOL, NROW, NLAY) ,  CC (NCOL, NROW, NLAY) , 

3  CV (NCOL, NROW, NLAY) ,  VBNM(4,20),  VBVL(4,20), 

4  SC2  (NCOL, NROW, NLAY)  , 

5  TOP (NCOL, NROW, NLAY) , BUFF (NCOL, NROW, NLAY) 

C 

COMMON  /FLWCOM/ LAYCON (80) 

C 

DIMENSION  TEXT (4) 

C 

DATA  TEXT (1) , TEXT (2) , TEXT (3) ,TEXT(4)  /'  STO','RAGE'/ 


C  - 

C 

d - INITIALIZE  BUDGET  ACCUMULATORS 

STOIN=0 . 

STOUT=0 . 

C 

C2 - IF  CELL-BY-CELL  FLOWS  ARE  NEEDED  THEN  SET  FLAG  IBD . 

IBD=0 

IF (ICBCFL.NE.O  .AND.  IGFDCB.GT.O)  IBD=1 
C 

C3 - IF  STEADY  STATE  THEN  SKIP  ALL  STORAGE  CALCULATIONS 

IF(ISS.NE.O)  GO  TO  305 
C 

C4 - IF  CELL-BY-CELL  FLOWS  ARE  NEEDED  (IBD  IS  SET)  CLEAR  BUFFER 

IF(IBD.EQ.O)  GO  TO  220 


DO  210  K=l, NLAY 
DO  210  1=1, NROW 
DO  210  J=1 , NCOL 
BUFF ( J, I , K) =0 . 

210  CONTINUE 
C 

C5 - RUN  THROUGH  EVERY  CELL  IN  THE  GRID 

220  KT=0 

DO  300  K=l, NLAY 
LC=LAYCON(K) 

IF  (LC . EQ . 3  .OR.  LC.EQ.2)  KT-KT  +  1 
DO  300  1=1, NROW 
DO  300  J=l, NCOL 
C 
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C6 - CALCULATE  FLOW  FROM  STORAGE  (VARIABLE  HEAD  CELLS  ONLY) 

IF ( IBOUND ( J, I , K) . LE . 0 )  GO  TO  300 
HSING=HNEW { J, I, K) 

C 

C6A - CHECK  LAYER  TYPE  TO  SEE  IF  ONE  STORAGE  FACTOR  OR  TWO 

IF (LC .NE . 3  .AND.  LC.NE.2)  GO  TO  285 
C 

C6B - TWO  STORAGE  CAPACITIES 


TP=TOP ( J, I , KT ) 

SYA=SC2 (J, I,KT) 

SCFA=SC1 (J, I,  K) 

SOLD=SYA 

IF(HOLD(J,I,K) .GT.TP)  SOLD=SCFA 
SNEW=SYA 

IF (HSING.GT .TP)  SNEW=SCFA 

STRG=SOLD* (HOLD ( J, I, K) -TP)  +  SNEW*TP  -  SNEW*HSING 
GO  TO  288 
C 


C6C - ONE  STORAGE  CAPACITY 

285  SC=SC1 (J, I,K) 

STRG=SC*HOLD ( J, I, K)  -  SC*HSING 
C 

C7 - STORE  CELL-BY-CELL  FLOW  IN  BUFFER  AND  ADD  TO  ACCUMULATORS 


288  IF(IBD.EQ.l)  BUFF (J, I,K) =STRG/DELT 
IF(STRG)  292,300,294 
292  STOUT=STOUT-STRG 
GO  TO  300 

294  STOIN=STOIN+STRG 
C 

300  CONTINUE 


C 

c8 - IF  IBD  flag  IS  SET  RECORD  THE  CONTENTS  OF  THE  BUFFER 

IF(IBD.EQ.l)  CALL  UBUDSV (KSTP, KPER, TEXT, 

1  IGFDCB, BUFF, NCOL, NROW, NLAY, IOUT) 

C 

09 - ADD  TOTAL  RATES  AND  VOLUMES  TO  VBVL  &  PUT  TITLES  IN  VBNM 


305  VBVL (1,MSUM) =VBVL ( 1 , MSUM) +STOIN 
VBVL (2, MSUM) =VBVL (2, MSUM) +STOUT 
VBVL ( 3 , MSUM) =STOIN/DELT 
VBVL ( 4 , MSUM) =STOUT/DELT 
VBNM ( 1 , MSUM) =TEXT ( 1 ) 

VBNM ( 2 , MSUM) =TEXT (2) 

VBNM ( 3 , MSUM) =TEXT ( 3 ) 

VBNM ( 4 , MSUM) =TEXT ( 4 ) 

MSUM=MSUM+ 1 


C 

CIO - CALCULATE  FLOW  FROM  CONSTANT  HEAD  NODES 

CALL  SGFD1F (VBNM, VBVL, MSUM, HNEW, IBOUND, CR, CC, CV, TOP, DELT, 

1  NCOL, NROW, NLAY, KSTP, KPER, IBD, IGFDCB, ICBCFL, BUFF, IOUT) 

C 

Cll - CALCULATE  AND  SAVE  FLOW  ACROSS  CELL  BOUNDARIES  IF  C-B-C 

Cll - FLOW  TERMS  ARE  REQUESTED. 

IF (IBD.NE.0)  CALL  SGFDlB (HNEW, IBOUND, CR, CC, CV, TOP , NCOL, NROW, NLAY, 
1  KSTP, KPER, IGFDCB, BUFF, IOUT) 

C 

C12 - RETURN 

RETURN 

END 
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List  of  Variables  for  Module  GFD1BD 


Variable 

Range 

Definition 

BUFF 

Global 

DIMENSION  (NCOL, NROW, NLAY)  ,  Buffer  used  to  accumulate 
information  before  printing  or  recording  it. 

CC 

Global 

DIMENSION  {NCOL, NROW, NLAY) ,  Conductance  in  the  column 
direction 

CR 

Global 

DIMENSION  {NCOL, NROW, NLAY) ,  Conductance  in  the  row 
direction . 

CV 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  vertical 
direction . 

DELT 

GLOBAL 

Length  of  the  current  time  step. 

HNEW 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Most  recent  estimate  of  head  in 
each  cell. 

HOLD 

Global 

DIMENSION  {NCOL, NROW, NLAY) ,  Head  at  the  start  of  the  current 
time  step. 

HSING 

Module 

Temporary  single  precision  equivalent  of  HNEW(J,I,K). 

I 

Module 

Index  for  rows . 

IBD 

Package 

Flag . 

=  0,  cell-by-cell  flow  terms  for  this  package  will  not 
be  recorded  this  time  step. 

=  1,  cell-by-cell  flow  terms  for  this  package  will  be 
recorded  this  time  step. 

I BOUND 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Status  of  each  cell: 

<  0,  constant-head  cell 
=  0,  no-flow  cell 
>  0,  variable-head  cell 

ICBCFL 

Global 

Flag . 

=  0,  cell-by-cell  flow  terms  will  not  be  recorded  or 
printed  this  time  step. 

not  0,  cell-by-cell  flow  terms  for  this  time  step  will 
be  either  printed  (flow  to  constant-head  cells) 
or  recorded  (all  GFD  terms)  depending  on  IGFDCB. 

IGFDCB 

Package 

Flag  and  a  unit  number. 

>  0,  unit  number  on  which  the  cell-by-cell  flow  terms 
will  be  recorded  when  ICBCFL  is  set . 

=  0,  cell-by-cell  flow  terms  for  GFD  will  not  be  printed 
or  recorded. 

<  0,  flow  from  individual  constant-head  cells  will  be 
printed  when  ICBCFL  is  set . 

I  OUT 

Global 

Primary  unit  number  for  all  printed  output . 

ISS 

Package 

Steady-state  flag. 

=  0,  simulation  is  transient 
not  0,  simulation  is  steady  state. 

J 

Module 

Index  for  columns . 

K 

Module 

Index  for  layers . 

KPER 

Global 

Stress  period  counter. 

KSTP 

Global 

Time  step  counter.  KSTP-1  at  the  start  of  each  stress 
period. 

KT 

Module 

Index  for  layer  within  TOP  and  (if  transient)  SC2  arrays. 

LAYCON 

Package 

DIMENSION  (80)  Layer-type  code: 

0  --  Layer  is  strictly  confined. 

1  --  Layer  is  strictly  unconfined. 

2  --  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  is  constant) . 

3  —  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  varies) . 

LC 

Module 

Temporary  equivalent  of  LAYCON (K) . 
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MSUM 

NCOL 

NLAY 

NROW 

SC 

SCI 

SC2 


SCFA 

SNEW 

SOLD 

STOIN 

STOUT 

STRG 

SYA 

TEXT 

TOP 


TP 

VBNM 

VBVL 


Global 

Global 

Global 

Global 

Module 

Package 

Package 


Module 

Module 

Module 

Module 

Module 

Module 

Module 

Module 

Package 


Module 

Global 

Global 


Index  for  budget  entries  in  VBVL  and  VBNM. 

Number  of  columns  in  the  grid. 

Number  of  layers  in  the  grid. 

Number  of  rows  in  the  grid. 

Temporary  equivalent  of  SC1(J,I,K). 

DIMENSION  (NCOL, NROW, NLAY) ,  Primary  storage  capacity. 
DIMENSION  (NCOL, NROW, NLAY) ,  Secondary  storage  capacity. 
Although  SC2  is  dimensioned  to  the  size  of  the  grid, 
space  exists  only  for  cells  that  can  convert  between 
confined  and  unconfined. 

Temporary  equivalent  of  SC1(J,I,K). 

Storage  capacity  at  the  end  of  the  time  step. 

Storage  capacity  at  the  start  of  the  time  step. 

Sum  of  decreases  in  storage  volume  for  individual  cells . 
Sum  of  increases  in  storage  volume  for  individual  cells. 
Volume  into  or  out  of  storage  for  an  individual  cell. 
Temporary  equivalent  of  SC2 ( J, I,K). 

Label  for  storage  term  in  volumetric  budget. 

DIMENSION  (NCOL, NROW, NLAY) ,  Elevation  of  the  aquifer  top. 
Although  TOP  is  dimensioned  to  the  size  of  the  grid, 
space  exists  only  for  cells  that  can  convert  between 
confined  and  unconfined. 

Temporary  equivalent  of  TOP(J,I,KT) . 

DIMENSION  (4,20),  Labels  for  terms  in  volumetric  budget. 
DIMENSION  (4,20),  Values  for  terms  in  volumetric  budget. 
For  term  "N",  the  values  are: 

(1,N),  Flow  rate  into  the  model  this  time  step. 

(2, N) ,  Flow  rate  out  of  the  model  this  time  step. 
(3,N) ,  Cumulative  volume  into  the  model. 

(4,N),  Cumulative  volume  out  of  the  model. 
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SGFD1N  Module 


SUBROUTINE  SGFD1N (HNEW, IBOUND, CR, CC, CV, CDTR, CDTC, 

1  NCOL , NROW, NLAY, IOUT ) 

C 

C - VERSION  1438  19SEP1989  SGFD1N 

C 

0  ****************************************************************** 

C  INITIALIZE  AND  CHECK  GFD  DATA 

c 

C  SPECIFICATIONS: 

C  - 

c 

DOUBLE  PRECISION  HNEW, HCNV 
C 

DIMENSION  HNEW {NCOL, NROW, NLAY) , IBOUND (NCOL, NROW, NLAY) , 

1  CDTR (NCOL, NROW, NLAY) , CR (NCOL, NROW, NLAY) , 

2  CC (NCOL, NROW, NLAY) , CV (NCOL, NROW, NLAY) , 

3  CDTC (NCOL, NROW, NLAY) 

C 

COMMON  / FLWCOM/ LAYCON (80) 


Cl - IF  IBOUND=0 ,  SET  CR=CC=CV=CDTR=CDTC=0 . 

KB=0 

DO  30  K=1 , NLAY 

IF (LAYCON (K) . EQ . 3  .OR.  LAYCON (K) . EQ . 1 )  KB=KB+1 
DO  30  1=1, NROW 
DO  30  J=1 , NCOL 

IF (IBOUND (J, I,K) .NE. 0)  GO  TO  30 
IF (K.NE .NLAY)  CV(J,I,K)=0. 

IF(K.NE.l)  CV(J, I,K-1)=0 . 

CC  ( J,  I ,  K)  =0  . 

IF(I.NE.l)  CC ( J, 1-1 , K) =0 . 

CR(J, I,K) =0 . 

IF(J.NE.l)  CR ( J-l, I, K) =0 . 

IF (LAYCON (K) . NE . 3  .AND.  LAYCON (K) . NE . 1 )  GO  TO  30 
CDTR ( J, I, KB) =0 . 

CDTC (J, I, KB) =0 . 

IF(J.NE.l)  CDTR (J-l, I, KB) =0 . 

IF(I.NE.l)  CDTC ( J, 1-1, KB) =0 . 

30  CONTINUE 
C 

C2 - CHECK  IF  ANY  ACTIVE  NODE  WILL  HAVE  ALL  ZERO  CONDUCTANCE 

C2 - IF  SO,  CONVERT  NODE  TO  NOFLOW 

HCNV=888 .88 
KB=0 

DO  70  K=l, NLAY 

IF (LAYCON (K) . EQ . 1  .OR.  LAYCON (K) . EQ . 3 )  GO  TO  55 

C2A - WHEN  LAYER  TYPE  0  OR  2,  CR,  CC,  OR  CV  MUST  BE  NONZERO 

DO  54  1=1, NROW 
DO  54  J=1 , NCOL 

IF (IBOUND (J, I,K) .EQ. 0)  GO  TO  54 
IF(J.EQ.l)  GO  TO  41 
IF (CR ( J-l , I , K) .NE.0.)  GO  TO  54 
41  IF (J.EQ.NCOL)  GO  TO  43 

IF(CR(J,I,K) .NE.0.)  GO  TO  54 
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43  IF(I.EQ.l)  GO  TO  45 

IF(CC(J,I-1,K) .NE.O.)  GO  TO  54 
45  IF(I.EQ.NROW)  GO  TO  47 

IF(CC(J, I,K) .NE.O.)  GO  TO  54 
47  IF (K.EQ.NLAY)  GO  TO  51 

IF(CV(J,I,K) .NE.O.)  GO  TO  54 

51  IF(K.EQ.l)  GO  TO  53 

IF (CV ( J, I, K-l) .NE.O.)  GO  TO  54 

53  IBOUND ( J, I, K) =0 
HNEW(J, I ,  K) =HCNV 
WRITE (IOUT, 52)  K,  I,J 

52  FORMAT (IX, 'NODE  ( LAYER, ROW, COL) r , 3 I 4 , 

1  '  ELIMINATED  BECAUSE  ALL  CONDUCTANCES  TO  NODE  ARE  O') 

54  CONTINUE 
GO  TO  70 

C 

C2B - WHEN  LAYER  TYPE  IS  1  OR  3,  CDTR,  CDTC,  OR  CV  MUST  BE  NONZERO 

55  KB=KB+1 

DO  69  1=1 , NROW 
DO  69  J=1 , NCOL 

IF(IBOUND(J,I,K) .EQ.O)  GO  TO  69 
IF (K.EQ.NLAY)  GO  TO  56 
IF <CV(J, I,K) .NE.O. )  GO  TO  69 

56  IF(K.EQ.l)  GO  TO  57 

IF (CV ( J, I, K-l) .NE.O.)  GO  TO  69 

57  IF(J.EQ.l)  GO  TO  59 

IF (CDTR ( J-l, I, KB) .NE.O.)  GO  TO  69 
59  IF ( J.EQ.NCOL)  GO  TO  61 

IF (CDTR (J, I, KB) .NE.O.)  GO  TO  69 
61  IF(I.EQ.l)  GO  TO  63 

IF (CDTC ( J, 1-1 , KB) .NE.O.)  GO  TO  69 
63  IF (I .EQ. NROW)  GO  TO  67 

IF(CDTC(J, I,KB) .NE.O.)  GO  TO  69 
67  IBOUND (J, I  ,K) =0 
HNEW ( J, I,K) =HCNV 
CC (J, I,K) =0 . 

WRITE (IOUT, 52)  K,  I,J 

69  CONTINUE 

70  CONTINUE 
C 

C6 - RETURN 

100  RETURN 
END 
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List  of  Variables  for  Module  SGFD1N 


Variable 

Range 

Definition 

CC 

Global 

DIMENSION  (NCOL, NROW, NLAY)  ,  Conductance  in  the  column 
direction . 

CDTC 

Package 

DIMENSION  (NCOL, NROW, NLAY)  ,  Conductance  divided  by 

thickness  in  the  column  direction.  Although  CDTC  is 
dimensioned  to  the  size  of  the  grid,  space  exists  only 
for  cells  for  which  saturated  thickness  is  calculated. 

CDTR 

Package 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  divided  by 
thickness  in  the  row  direction.  Although  CDTR  is 
dimensioned  to  the  size  of  the  grid,  space  exists  only 
for  cells  for  which  saturated  thickness  is  calculated. 

CR 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  row 
direction . 

CV 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  vertical 
direction . 

HCNV 

Module 

Value  of  head  used  to  indicate  that  a  cell  has  been 

converted  to  no  flow  because  all  conductances  to  that 
cell  are  0 . 

HNEW 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Most  recent  estimate  of  head  in 
each  cell. 

I 

Module 

Index  for  rows . 

I BOUND 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Status  of  each  cell: 

<  0,  constant-head  cell 
=  0,  no-flow  cell 
>  0,  variable-head  cell 

I  OUT 

Global 

Primary  unit  number  for  all  printed  output. 

J 

Module 

Index  for  columns . 

K 

Module 

Index  for  layers . 

KB 

Module 

Index  for  layer  within  CDTR  and  CDTC  arrays. 

LAYCON 

Package 

DIMENSION  (80)  Layer-type  code: 

0  —  Layer  is  strictly  confined. 

1  —  Layer  is  strictly  unconfined. 

2  —  Layer  ‘is  convertible  between  confined  and  unconfined 

(transmissivity  is  constant) . 

3  —  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  varies) . 

NCOL 

Global 

Number  of  columns  in  the  grid. 

NLAY 

Global 

Number  of  layers  in  the  grid. 

NROW 

Global 

Number  of  rows  in  the  grid. 
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SGFD1H  Module 


SUBROUTINE  SGFD1H (HNEW, IBOUND, CR, CC, CV, CDTR, CDTC, BOT, 

1  TOP , K, KB, KT, KITER, KSTP , KPER, NCOL, NROW,  NLAY, IOUT) 

C 

C - VERSION  1318  19SEP1989  SGFD1H 

C 

Q  A***************************************************************** 

C  COMPUTE  CONDUCTANCE  FROM  SATURATED  THICKNESS  AND  CONDUCTANCE 

C  DIVIDED  BY  THICKNESS 

q  ****************************************************************** 

C 

C  SPECIFICATIONS: 

C  - 

DOUBLE  PRECISION  HNEW 
C 

DIMENSION  HNEW (NCOL, NROW, NLAY) , IBOUND (NCOL, NROW, NLAY) , 

1  CR (NCOL, NROW, NLAY) ,  CC (NCOL, NROW, NLAY) ,  CV (NCOL, NROW, NLAY) , 

2  CDTR (NCOL, NROW, NLAY) ,  CDTC (NCOL, NROW, NLAY) , 

3  BOT (NCOL, NROW, NLAY) , TOP (NCOL, NROW, NLAY) 

C 

COMMON  /FLWCOM/ LAYCON (80) 

C  - 

C 

Cl - CALCULATE  SATURATED  THICKNESS  AT  EACH  ACTIVE  NODE  AND  STORE 

C1 - TEMPORARILY  IN  CC 

DO  200  1=1, NROW 
DO  200  J=1 , NCOL 
C 

CIA - IF  CELL  IS  INACTIVE  THEN  SET  THICKNESS  =  0. 

IF (IBOUND (J, I,K) .NE. 0)  GO  TO  10 
CC (J, I , K) =0 . 

GO  TO  200 
C 

C1B - CALCULATE  SATURATED  THICKNESS. 

10  HD=HNEW ( J, I, K) 

IF (LAYCON (K) .EQ. 1)  GO  TO  50 
IF (HD.GT.TOP (J, I, KT) )  HD=TOP (J, I,KT) 

50  THCK=HD-BOT ( J, I, KB) 

IF (THCK.LE.O . )  GO  TO  100 
C 

C1C - IF  SATURATED  THICKNESS>0  THEN  SAVE  IT  IN  CC 

CC(J, I,K)=THCK 
GO  TO  200 
C 

C1D - WHEN  SATURATED  THICKNESS  <  0,  PRINT  A  MESSAGE  AND  SET 

C1D - IBOUND,  THICKNESS,  AND  VERTICAL  CONDUCTANCE  =0 

100  WRITE (IOUT, 150)  K, I , J, KITER, KSTP , KPER 

150  FORMAT (1H0, 10 (' *' ) , 'NODE' , 314, '  ( LAYER, ROW, COL)  WENT  DRY' 

1  , '  AT  ITERATION  =  ' , 13, '  TIME  STEP  =' , 13 

2  ,'  STRESS  PERIOD  =',I3) 

HNEW ( J, I , K) =1 .E30 

CC  (J,  I ,  K)  =0  . 

IBOUND ( J, I, K) =0 
IF (K.LT.NLAY)  CV(J,I,K)=0. 

IF(K.GT.l)  CV ( J, I, K-l) =0 . 

200  CONTINUE 
C 
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C 2 - COMPUTE  HORIZONTAL  BRANCH  CONDUCTANCE  FROM  SATURATED 

C2 - THICKNESS  AT  NODES  AND  CDT  BETWEEN  NODES. 

C2 - FOR  EACH  CELL  CALCULATE  BRANCH  CONDUCTANCE  FROM  THAT  CELL 

C2 - TO  THE  ONE  ON  THE  RIGHT  AND  THE  ONE  IN  FRONT. 


DO  400  1=1 , NROW 
DO  400  J=l, NCOL 
B1=CC(J, I,K) 

C 

C2A - IF  B1=0  THEN  SET  CONDUCTANCE  EQUAL  TO  0 .  GO  ON  TO  NEXT  CELL. 

IF (B1 . NE . 0 . )  GO  TO  310 
CR(J, I,K)=0 . 

GO  TO  400 


C 

C2B - IF  THIS  IS  NOT  THE  LAST  COLUMN (RIGHTMOST)  THEN  CALCULATE 

C2B - BRANCH  CONDUCTANCE  IN  THE  ROW  DIRECTION  (CR)  TO  THE  RIGHT. 


310  IF(J.EQ.NCOL)  GO  TO  330 
62=00(0+1, I ,  K) 

RATIO=B2/Bl 

IF (RATIO. NE.0.)  GO  TO  320 
CR(J,I,K)=0. 

GO  TO  330 
C 

320  THCK= (B1+B2) * .5 

IF (RATIO . GE .1.25  .OR.  RATIO . LE . 0 . 8 )  THCK= (B2-B1) /ALOG (RATIO) 
CR(J, I, K)=THCK*CDTR(J, I, KB) 


C 

C2C - IF  THIS  IS  NOT  THE  LAST  ROW (FRONTMOST)  THEN  CALCULATE 

C2C - BRANCH  CONDUCTANCE  IN  THE  COLUMN  DIRECTION  (CC)  TO  THE  FRONT. 


330  IF (I. EQ. NROW)  GO  TO  400 
B2=CC ( Jf 1+1, K) 

RATIO=B2/Bl 

IF (RATIO. NE.0 . )  GO  TO  340 
CC(J,I,K)=0. 

GO  TO  400 
C 

340  THCK= (B1+B2 ) * . 5 

IF (RATIO. GE. 1.25  .OR.  RATIO . LT . 0 . 8)  THCK= (B2-B1 ) /ALOG (RATIO) 
CC  ( J,  I,  K) =THCK*CDTC ( J, I, KB) 

400  CONTINUE 
C 

C3 - RETURN 

RETURN 

END 
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List  of  Variable  for  Module  SGFD1H 


Variable 


Bl 

B2 

BOT 


CC 


CDTC 


CDTR 


CR 

CV 

HD 

HNEW 

I 

I BOUND 


IOUT 

J 

K 

KB 

KITER 

KPER 

KSTP 

KT 

LAYCON 


NCOL 

NLAY 

NROW 

RATIO 

THCK 

TOP 


Range 


Definition 


Module 

Module 

Package 


Global 


Package 


Package 


Global 

Global 

Module 

Global 

Module 

Global 


Global 

Module 

Module 

Module 

Global 

Global 

Global 

Module 

Package 


Global 

Global 

Global 

Module 

Module 

Package 


Saturated  thickness  at  node  (J, I,K)  . 

Saturated  thickness  at  node  (J+1,I,K)  or  node  (J, I+1,K). 

DIMENSION  (NCOL, NROW, NLAY) ,  Elevation  of  the  aquifer 

bottom.  Although  BOT  is  dimensioned  to  the  size  of  the 
grid,  space  exists  only  for  cells  for  which  saturated 
thickness  is  calculated. 


DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  column 

direction.  This  array  is  also  used  to  temporarily  hold 
saturated  thickness. 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  divided  by 

thickness  in  the  column  direction.  Although  CDTC  is 
dimensioned  to  the  size  of  the  grid,  space  exists  only 
for  cells  for  which  saturated  thickness  is  calculated. 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  divided  by 
thickness  in  the  row  direction.  Although  CDTR  is 
dimensioned  to  the  size  of  the  grid,  space  exists  only 
for  cells  for  which  saturated  thickness  is  calculated. 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  row 
direction . 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  vertical 
direction . 

Single  precision  equivalent  of  HNEW(J,I,K). 

DIMENSION  (NCOL, NROW, NLAY) ,  Most  recent  estimate  of  head  in 
each  cell. 

Index  for  rows. 


DIMENSION  (NCOL, NROW, NLAY) ,  Status  of  each  cell: 

<  0,  constant-head  cell 
=  0,  no-flow  cell 
>  0,  variable-head  cell 
Primary  unit  number  for  all  printed  output. 

Index  for  columns. 

Index  for  layers. 

Index  for  layer  within  BOT,  CDTR,  and  CDTC  arrays. 

Iteration  counter.  KITER=1  at  the  start  of  each  time  step. 
Stress  period  counter. 

Time  step  counter.  KSTP=1  at  the  start  of  each  stress 
period. 

Index  for  layer  within  TOP  array. 

DIMENSION  (80)  Layer-type  code: 

0  --  Layer  is  strictly  confined. 

1  —  Layer  is  strictly  unconfined. 

2  --  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  is  constant) . 

3  --  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  varies) . 

Number  of  columns  in  the  grid. 

Number  of  layers  in  the  grid. 

Number  of  rows  in  the  grid. 

B2/B1 . 

Saturated  thickness. 


DIMENSION  (NCOL, NROW, NLAY) ,  Elevation  of  the  aquifer  top. 
Although  TOP  is  dimensioned  to  the  size  of  the  grid, 
space  exists  only  for  cells  that  can  convert  between 
confined  and  unconfined. 
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SGFDlB  Module 


SUBROUTINE  SGFDlB (HNEW, I BOUND, CR, CC, CV, TOP , NCOL, NROW, NLAY, 

1  KSTP, KPER, IGFDCB, BUFF, IOUT) 

C 

c - VERSION  1328  19SEP1989  SGFDlB 

C 

0  ****************************************************************** 

C  COMPUTE  FLOW  ACROSS  EACH  CELL  WALL 

0  ****************************************************************** 

C 

C  SPECIFICATIONS: 

C  - 

CHARACTER* 4  TEXT 
DOUBLE  PRECISION  HNEW, HD 
C 

DIMENSION  HNEW (NCOL, NROW, NLAY) ,  IBOUND (NCOL, NROW, NLAY) , 

1  CR (NCOL, NROW, NLAY) ,  CC (NCOL, NROW, NLAY) , 

2  CV (NCOL, NROW, NLAY) ,  TOP (NCOL, NROW, NLAY) , 

3  BUFF (NCOL, NROW,  NLAY) 

C 

COMMON  / F  LWCOM / LAYCON (80) 

C 

DIMENSION  TEXT (12) 

C 

DATA  TEXT ( 1 ) , TEXT (2) , TEXT (3) , TEXT (4) , TEXT (5) , TEXT (6) , TEXT (7) , 

1  TEXT ( 8 ) , TEXT ( 9 ) ,TEXT(10) ,TEXT(11) ,TEXT(12) 

2  /'FLOW','  RIG' , ' HT  F','ACE  ', 

2  'FLOW','  FRO', 'NT  F','ACE  ','FLOW','  LOW' , ' ER  F' , ' ACE  '/ 

0  - 

C 

NCMl=NCOL- 1 

IF(NCMl.LT.l)  GO  TO  405 
C 

01 - CLEAR  THE  BUFFER 

DO  310  K=1 , NLAY 
DO  310  1=1, NROW 
DO  310  J=1 , NCOL 
BUFF ( J, I, K) =0 . 

310  CONTINUE 
C 

C2 - FOR  EACH  CELL  CALCULATE  FLOW  THRU  RIGHT  FACE  &  STORE  IN  BUFFER 

DO  400  K=l, NLAY 
DO  400  1=1, NROW 
DO  400  J=1 , NCM1 

IF ( ( IBOUND ( J, I , K) . LE . 0 )  .AND.  ( IBOUND (J+l , I , K) . LE . 0 ) )  GO  TO  400 
HDIFF=HNEW ( J, I,K) -HNEW (J+l, I,K) 

BUFF (J, I , K) =HDIFF*CR ( J, I,K) 

400  CONTINUE 
C 

C3 - RECORD  CONTENTS  OF  BUFFER 

CALL  UBUDSV (KSTP, KPER, TEXT (1) , IGFDCB, BUFF, NCOL, NROW, NLAY, IOUT) 

C 

04 - CLEAR  THE  BUFFER 

405  NRM 1 =NROW - 1 

IF (NRM1 .LT. 1)  GO  TO  505 
DO  410  K=l, NLAY 
DO  410  1=1, NROW 
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DO  410  J=l,NCOL 
BUFF (J, I,K)=0 . 

410  CONTINUE 
C 

C5 - FOR  EACH  CELL  CALCULATE  FLOW  THRU  FRONT  FACE  &  STORE  IN  BUFFER 

DO  500  K=l, NLAY 
DO  500  I=1/NRM1 
DO  500  J=l,NCOL 

IF ( ( I BOUND { J,  I ,  K) . LE . 0 )  .AND.  ( I BOUND (J, I+1,K) .LE.0) )  GO  TO  500 
HDIFF=HNEW(J, I , K) -HNEW(J, I+1,K) 

BUFF ( J, I, K) =HDIFF*CC ( J, I, K) 

500  CONTINUE 
C 

C6 - RECORD  CONTENTS  OF  BUFFER. 

CALL  UBUDSV (KSTP, KPER, TEXT (5) , IGFDCB, BUFF, NCOL, NROW, NLAY, IOUT) 
505  NLM1=NLAY- 1 

IF(NLMl.LT.l)  GO  TO  1000 
C 

C7 - CLEAR  THE  BUFFER 

DO  510  K=1 , NLAY 
DO  510  1=1, NROW 
DO  510  J=1 , NCOL 
BUFF {J, I,K)=0 . 

510  CONTINUE 
C 

C8-- - FOR  EACH  CELL  CALCULATE  FLOW  THRU  LOWER  FACE  &  STORE  IN  BUFFER 

KT=0 

DO  600  K=1 , NLM1 

IF(LAYCON(K) . EQ . 3  .OR.  LAYCON (K) . EQ . 2 )  KT=KT+1 
DO  600  1=1 , NROW 
DO  600  J=l, NCOL 

IF ( (IBOUND (J, I,K) .LE.0)  .AND.  (IBOUND ( J, I, K+l) . LE . 0) )  GO  TO  600 
HD=HNEW ( J, I, K+l) 

IF (LAYCON (K+l) .NE.3  .AND.  LAYCON (K+l ). NE . 2 )  GO  TO  580 
TMP=HD 

IF (TMP . LT . TOP (J, I , KT+1 ) )  HD=TOP (J, I, KT+1) 

580  HDIFF=HNEW(J, I,K) -HD 

BUFF (J, I , K) =HDIFF*CV ( J, I,K) 

600  CONTINUE 
C 

C9 - RECORD  CONTENTS  OF  BUFFER. 

CALL  UBUDSV (KSTP , KPER, TEXT ( 9 ) , IGFDCB, BUFF, NCOL, NROW, NLAY, IOUT) 

C 

CIO - RETURN 

1000  RETURN 
END 
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List  of  Variable  for  Module  SGFD1B 


Variable 

Range 

Definition 

BUFF 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Buffer  used  to  accumulate 
information  before  printing  or  recording  it. 

CC 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  column 
direction . 

CR 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  row 
direction . 

CV 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  vertical 
direction . 

HD 

Module 

Temporary  equivalent  of  HNEW ( J, I, K+l)  or  TOP ( J, I, K+l) . 

HDIFF 

Module 

Head  difference  between  two  adjacent  nodes. 

HNEW 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Most  recent  estimate  of  head  in 
each  cell. 

I 

Module 

Index  for  rows . 

I BOUND 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Status  of  each  cell: 

<  0,  constant-head  cell 
=  0,  no-flow  cell 
>  0,  variable-head  cell 

IGFDCB 

Package 

Flag  and  a  unit  number. 

>  0,  unit  number  on  which  the  cell-by-cell  flow  terms 
will  be  recorded  when  ICBCFL  is  set . 

=  0,  cell-by-cell  flow  terms  for  GFD  will  not  be  printed 
or  recorded. 

<  0,  flow  from  individual  constant-head  cells  will  be 
printed  when  ICBCFL  is  set. 

IOUT 

Global 

Primary  unit  number  for  all  printed  output. 

J 

Module 

Index  for  columns. 

K 

Module 

Index  for  layers. 

KPER 

Global 

Stress  period  counter. 

KSTP 

Global 

Time  step  counter.  KSTP=1  at  the  start  of  each  stress 
period. 

KT 

Module 

Index  for  layer  within  TOP  array. 

LAYCON 

Package 

DIMENSION  (80)  Layer-type  code: 

0  —  Layer  is  strictly  confined. 

1  —  Layer  is  strictly  unconfined. 

2  —  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  is  constant) . 

3  —  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  varies) . 

NCMl 

Module 

NCOL-1 . 

NCOL 

Global 

Number  of  columns  in  the  grid. 

NLAY 

Global 

Number  of  layers  in  the  grid. 

NLM1 

Module 

NLAY-1 

NRM1 

Module 

NROW-1 . 

NROW 

Global 

Number  of  rows  in  the  grid. 

TEXT 

Module 

Labels  for  flow  across  cell  walls  in  cell-by-cell  budget 
file . 

TMP 

Module 

Single  precision  equivalent  of  HD. 

TOP 

Package 

DIMENSION  (NCOL, NROW, NLAY) ,  Elevation  of  the  aquifer  top. 
Although  TOP  is  dimensioned  to  the  size  of  the  grid, 
space  exists  only  for  cells  that  can  convert  between 
confined  and  unconfined. 
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SGFD1F  Module 


SUBROUTINE  SGFDlF (VBNM, VBVL,MSUM, HNEW, IBOUND, CR, CC, CV, 

1  TOP , DELT, NCOL, NROW, NLAY, KSTP, KPER, IBD, IGFDCB, ICBCFL, 

2  BUFF, IOUT) 

C - VERSION  1429  19SEP1989  SGFDlF 

C 

C  ****************************************************************** 

C  COMPUTE  FLOW  FROM  CONSTANT  HEAD  NODES 

C  ******************************************************************** 

c 

C  SPECIFICATIONS: 

C  - 

CHARACTER* 4  VBNM, TEXT 
DOUBLE  PRECISION  HNEW, HD 
C 

DIMENSION  HNEW (NCOL, NROW, NLAY) ,  IBOUND (NCOL, NROW, NLAY) , 

1  CR (NCOL, NROW, NLAY) ,  CC (NCOL, NROW, NLAY) , 

2  CV (NCOL, NROW, NLAY) ,  VBNM (4, 20),  VBVL(4,20), 

3  TOP (NCOL, NROW, NLAY) , BUFF (NCOL, NROW, NLAY) 

C 

COMMON  /FLWCOM/LAYCON (80) 

C 

DIMENSION  TEXT (4) 

C 

DATA  TEXT (1) , TEXT (2) , TEXT (3) , TEXT (4)  /'  C' , ' ONST' , ' ANT  ','HEAD'/ 
C  - 

c 


Cl - CLEAR  BUDGET  ACCUMULATORS 

CHIN=0 . 

CHOUT=0 . 

C 

C2 - CLEAR  BUFFER  IF  CELL-BY-CELL  FLOW  TERM  FLAG (IBD)  IS  SET 

IF(IBD.EQ.O)  GO  TO  8 


DO  5  K=1 , NLAY 
DO  5  1=1, NROW 
DO  5  J=1 , NCOL 
BUFF (J, I, K) =0 . 
5  CONTINUE 


C 

C3 - FOR  EACH  CELL  IF  IT  IS  CONSTANT  HEAD  COMPUTE  FLOW  ACROSS  6 

C3 - FACES  . 


8  KT=0 

DO  200  K=1 , NLAY 
LC=LAYCON (K) 

IF (LC . EQ . 3  .OR.  LC.EQ.2)  KT=KT+1 
DO  200  1=1, NROW 
DO  200  J=1 , NCOL 


C 

C4 - IF  CELL  IS  NOT  CONSTANT  HEAD  SKIP  IT  &  GO  ON  TO  NEXT  CELL. 

IF  (IBOUND (J, I,K) .GE.O)GO  TO  200 
C 

C5 - CLEAR  FIELDS  FOR  SIX  FLOW  RATES. 


X1=0  . 
X2=0  . 
X3=0  . 
X4=0  . 
X5=0  . 
X6=0  . 
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C6 - FOR  EACH  FACE  OF  THE  CELL  CALCULATE  FLOW  THROUGH  THAT  FACE 

C6 - OUT  OF  THE  CONSTANT  HEAD  CELL  AND  INTO  THE  FLOW  DOMAIN. 

C6 - COMMENTS  7-11  APPEAR  ONLY  IN  THE  SECTION  HEADED  BY  COMMENT  6A 

C6 - BUT  THEY  APPLY  IN  A  SIMILAR  MANNER  TO  THE  SECTIONS  HEADED 

C6 - BY  COMMENTS  6B-6F. 

C 

C6A - CALCULATE  FLOW  THROUGH  THE  LEFT  FACE 

C 

C7 - IF  THERE  IS  NOT  A  VARIABLE  HEAD  CELL  ON  THE  OTHER  SIDE  OF  THIS 

C7 - FACE  THEN  GO  ON  TO  THE  NEXT  FACE. 

IF(J.EQ.l)  GO  TO  30 

IF (IBOUND (J-l, I,K) .LE. 0) GO  TO  30 

HDIFF^HNEW { J, I,K) -HNEW(J-1, I,K) 

C 

C8 - CALCULATE  FLOW  THROUGH  THIS  FACE  INTO  THE  ADJACENT  CELL. 

Xl=HDIFF*CR( J-l, I, K) 

C 

C9 - TEST  TO  SEE  IF  FLOW  IS  POSITIVE  OR  NEGATIVE 

IF  (XI)  10,30,20 
C 

CIO - IF  NEGATIVE  ADD  TO  CHOUT (FLOW  OUT  OF  DOMAIN  TO  CONSTANT  HEAD) . 

10  CHOUT=CHOUT-Xl 
GO  TO  30 
C 

Cll - IF  POSITIVE  ADD  TO  CHIN (FLOW  INTO  DOMAIN  FROM  CONSTANT  HEAD) . 

20  CHIN=CHIN+X1 
C 

C6B - CALCULATE  FLOW  THROUGH  THE  RIGHT  FACE 

30  IF (J.EQ.NCOL)  GO  TO  60 

IF (IBOUND (J+l, I,K) .LE. 0)  GO  TO  60 
HDIFF=HNEW(J, I, K) -HNEW ( J+l, I,K) 

X2=HDIFF*CR ( J, I , K) 

IF (X2) 40, 60,50 
40  CHOUT=CHOUT-X2 
GO  TO  60 
50  CHIN=CHIN+X2 
C 

C6C - CALCULATE  FLOW  THROUGH  THE  BACK  FACE. 

60  IF  (I .EQ. 1)  GO  TO  90 

IF  (IBOUND (J, I-1,K) .LE. 0)  GO  TO  90 
HDIFF=HNEW ( J, I,K) -HNEW(J, I-1,K) 

X3=HDIFF*CC ( J, 1-1 , K) 

IF (X3)  70,90,80 
70  CHOUT=CHOUT-X3 
GO  TO  90 
80  CHIN=CHIN+X3 
C 

C6D - CALCULATE  FLOW  THROUGH  THE  FRONT  FACE. 

90  IF (I .EQ.NROW)  GO  TO  120 

IF (IBOUND (J, I+1,K) .LE.0)  GO  TO  120 
HDIFF=HNEW ( J, I,K) -HNEW(J, I+1,K) 

X4=HDIFF*CC (J, I, K) 

IF  (X4 )  100,120,110 
100  CHOUT=CHOUT-X4 
GO  TO  120 
110  CHIN=CHIN+X4 
C 
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C6E - CALCULATE  FLOW  THROUGH  THE  UPPER  FACE 

120  IF(K.EQ.l)  GO  TO  150 

IF  (IBOUND (J, I, K-l) .LE.O)  GO  TO  150 
HD=HNEW(J, I,K) 

IF (LC .NE . 3  .AND.  LC.NE.2)  GO  TO  122 
TMP=HD 

IF (TMP . LT . TOP (J, I,KT) )  HD=TOP (J, I,KT) 

122  HDIFF=HD-HNEW(Jr I, K-l) 

X5=HDIFF*CV(J, I, K-l) 

IF (X5)  130,150,140 
130  CHOUT=CHOUT-X5 
GO  TO  150 
140  CHIN=CHIN+X5 
C 

C6F - CALCULATE  FLOW  THROUGH  THE  LOWER  FACE. 

150  IF (K.EQ.NLAY)  GO  TO  180 

IF (IBOUND (J, I, K+l) .LE.O)  GO  TO  180 
HD=HNEW ( J, I, K+l) 

IF ( LAYCON ( K+ 1 ) . NE . 3  .AND.  LAYCON (K+l) . NE . 2 )  GO  TO  152 
TMP=HD 

IF (TMP .LT. TOP (J, I , KT+1 ) )  HD=TOP (J, I,KT+1) 

152  HDIFF=HNEW ( J, I, K) -HD 
X6=HDIFF*CV ( J, I, K) 

IF (X6)  160,180,170 
160  CHOUT=CHOUT-X6 
GO  TO  180 
170  CHIN=CHIN+X6 
C 


C12 - SUM  UP  FLOWS  THROUGH  SIX  SIDES  OF  CONSTANT  HEAD  CELL. 

180  RATE=X1+X2+X3+X4+X5+X6 

C 

C13 - PRINT  THE  INDIVIDUAL  RATES  IF  REQUESTED ( IGFDCB<0 ) . 

IF ( IGFDCB . LT . 0 .AND . ICBCFL . NE . 0 )  WRITE ( IOUT, 900 )  (TEXT (N) ,N=1, 4) , 
1  KPER, KSTP, K, I, J, RATE 

900  FORMAT (1H0, 4A4, '  PERIOD' , 13, '  STEP', 13, r  LAYER', 13, 

1  '  ROW', 14,'  COL', 14,'  RATE  ',G15.7) 

C 


Cl 4 - IF  CELL-BY-CELL  FLAG  SET  STORE  SUM  OF  FLOWS  FOR  CELL  IN  BUFFER 

IF(IBD.EQ.l)  BUFF ( J, I, K) =RATE 
C 

200  CONTINUE 
C 

C15 - IF  CELL-BY-CELL  FLAG  SET  THEN  RECORD  CONTENTS  OF  BUFFER 

IF(IBD.EQ.l)  CALL  UBUDSV (KSTP , KPER, TEXT ( 1 ) , 

1  IGFDCB, BUFF, NCOL, NROW, NLAY, IOUT) 

C 

C 

Cl 6 - SAVE  TOTAL  CONSTANT  HEAD  FLOWS  AND  VOLUMES  IN  VBVL  TABLE 

C16 - FOR  INCLUSION  IN  BUDGET.  PUT  LABELS  IN  VBNM  TABLE. 

VBVL ( 1 , MSUM) =VBVL ( 1 , MSUM) +CHIN*DELT 
VBVL ( 2 , MSUM) =VBVL ( 2 , MSUM) +CHOUT *DELT 
VBVL ( 3 , MSUM) =CHIN 
VBVL ( 4 , MSUM) =CHOUT 
C 

C  - SETUP  VOLUMETRIC  BUDGET  NAMES 

VBNM ( 1 , MSUM) =TEXT ( 1 ) 

VBNM ( 2 , MSUM) =TEXT ( 2 ) 

VBNM ( 3 , MSUM) =TEXT ( 3 ) 

VBNM ( 4 , MSUM) =TEXT ( 4 ) 
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o  o 


c 

MSUM=MSUM+1 


Cl  7 - RETURN 

RETURN 

END 
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List  of  Variables  for  Module  SGFDlF 


Variable 

Range 

Definition 

BUFF 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Buffer  used  to  accumulate 
information  before  printing  or  recording  it. 

CC 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  column 
direction . 

CHIN 

Module 

Accumulator  for  flow  into  the  model  from  constant-head 
cells . 

CHOUT 

Module 

Accumulator  for  flow  out  of  the  model  to  constant-head 
cells . 

CR 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  row 
direction . 

CV 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Conductance  in  the  vertical 
direction . 

DELT 

GLOBAL 

Length  of  the  current  time  step. 

HD 

Module 

Temporary  equivalent  for  HNEW ( J, I, K+l)  or  TOP ( J, I , K+l ) . 

HD  IFF 

Module 

Head  difference  between  two  adjacent  nodes. 

HNEW 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Most  recent  estimate  of  head  in 
each  cell. 

I 

Module 

Index  for  rows . 

IBD 

Package 

Flag. 

=  0,  cell-by-cell  flow  terms  for  this  package  will  not 
be  recorded  this  time  step. 

=  1,  cell-by-cell  flow  terms  for  this  package  will  be 
recorded  this  time  step. 

I BOUND 

Global 

DIMENSION  (NCOL, NROW, NLAY) ,  Status  of  each  cell: 

<  0,  constant-head  cell 
=  0,  no-flow  cell 
>  0,  variable-head  cell 

ICBCFL 

Global 

Flag. 

=  0,  cell-by-cell  flow  terms  will  not  be  recorded  or 
printed  this  time  step. 

not  0,  cell-by-cell  flow  terms  for  this  time  step  will 
be  either  printed  (flow  to  constant-head  cells) 
or  recorded  (all  GFD  terms)  depending  on  IGFDCB. 

IGFDCB 

Package 

Flag  and  a  unit  number. 

>  0,  unit  number  on  which  the  cell-by-cell  flow  terms 
will  be  recorded  when  ICBCFL  is  set. 

=  0,  cell-by-cell  flow  terms  for  GFD  will  not  be  printed 
or  recorded. 

<  0,  flow  from  individual  constant-head  cells  will  be 
printed  when  ICBCFL  is  set . 

IOUT 

Global 

Primary  unit  number  for  all  printed  output . 

J 

Module 

Index  for  columns . 

K 

Module 

Index  for  layers . 

KPER 

Global 

Stress  period  counter. 

KSTP 

Global 

Time  step  counter.  KSTP=1  at  the  start  of  each  stress 
period. 

KT 

Module 

Index  for  layer  within  TOP  array. 

LAYCON 

Package 

DIMENSION  (80)  Layer-type  code: 

0  —  Layer  is  strictly  confined. 

1  —  Layer  is  strictly  unconfined. 

2  —  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  is  constant) . 

3  --  Layer  is  convertible  between  confined  and  unconfined 

(transmissivity  varies) . 

LC 

Module 

Temporary  equivalent  of  LAYCON (K) . 
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MSUM  Global  Index  for  budget  entries  in  VBVL  and  VBNM. 

N  Module  Index  for  TEXT  when  printing  cell-by-cell  constant-head 

flow. 

NCOL  Global  Number  of  columns  in  the  grid. 

NLAY  Global  Number  of  layers  in  the  grid. 

NROW  Global  Number  of  rows  in  the  grid. 

RATE  Module  Flow  from  a  constant-head  cell  into  the  model. 

TEXT  Module  Label  for  constant-head  flow  in  volumetric  budget. 

TMP  Module  Single  precision  equivalent  of  HD. 

TOP  Package  DIMENSION  (NCOL, NROW, NLAY) ,  Elevation  of  the  aquifer  top. 

Although  TOP  is  dimensioned  to  the  size  of  the  grid, 
space  exists  only  for  cells  that  can  convert  between 
confined  and  unconfined. 

VBNM  Global  DIMENSION  (4,20),  Labels  for  terms  in  volumetric  budget. 

VBVL  Global  DIMENSION  (4,20),  Values  for  terms  in  volumetric  budget. 

For  term  "N",  the  values  are: 

( 1 , N) ,  Flow  rate  into  the  model  this  time  step. 
(2,N),  Flow  rate  out  of  the  model  this  time  step. 
(3,N),  Cumulative  volume  into  the  model. 

(4,N),  Cumulative  volume  out  of  the  model. 

XI  Module  Flow  through  the  left  face. 

X2  Module  Flow  through  the  right  face. 

X3  Module  Flow  through  the  back  face. 

X4  Module  Flow  through  the  front  face. 

X5  Module  Flow  through  the  upper  face . 

X6  Module  Flow  through  the  lower  face. 
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